From 50d90d7aa42f60324a7fb9ba9dc45d4e2cac9d8e Mon Sep 17 00:00:00 2001 From: Chad Scherrer Date: Tue, 24 May 2022 09:46:17 -0700 Subject: [PATCH] Dev (#193) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Dirichlet(k::Integer, α) = Dirichlet(Fill(α, k)) * export TransformVariables as TV * drop redundant import * 0.0 => zero(Float64) * drop outdated Dists.logpdf * update StudentT * drop redundant import * update Uniform * bump MeasureBase version * reworking beta * small update to StudentT * basemeasure for discrete Distributions * using LogExpFunctions => import LogExpFunctions * quoteof(::Chain) * prettyprinting and chain-mucking * Some refactoring for Markov chains * import MeasureBase: ≪ * version bound for PrettyPrinting * copy(rng) might change its type (e.g. GLOBAL_RNG) * tests pass * cleaning up * more cleanup * big update * get tests passing * formatting * oops typo * move affine to MeasureTheory * updating * Val => StaticSymbol * more fixes * fix fix fix * more logdesnity => logdensity_def * more logdesnity fixes * debugging * formatting * bugfixes * working on tests * updates * working on tests * tests passing! * refactor * working on tests * drop static weight for now * fix sampling from ProductMeasure{<:Base.Generator} * tests passing!! * more stuff * constructor => constructorof * constructor =? construtorof * updates * working on tests * fix Dirichlet * update Bernoulli * working on tests * bugfixes for RealizedSamples * tests passing!! * tighten down inference * as(::PowerMeasure) * drop type-level stuff * using InverseFunctions.jl * update license * affero * copyright * update CI to 1.6 * xform => TV.as * oops missed a conflict * fix merge corruption * typo * fix license * Update README.md * merge * enumerate instead of zip * bugfix * inline rand * drop `static` from `insupport` results * update proxies * Move ConditionalMeasure to MeasureBase * IfElse.ifelse(p::Bernoulli, t, f) * IfElseMeasure * update some base measures * test broken :( * fix some redundancies * instance_type => Core.Typeof * update testvalue for Bernoulli and Binomial * un-break broken test (now passing) * Fall-back `For` method for when inference fails * drop extra spaces * more whitespace * bump MeasureBase dependency version * add newline * tidy up * ifelse tests * OEF newline * avoid type piracy * add Julia 1.7 to CI * make Julia 1.6 happy * approx instead of == * Require at least Julia 1.6 * Try Sebastian's idea test_measures ::Any[] * Another Any[] * Drop Likelihood test * drop 1.7 CI (seems buggy?) * bump version * export likelihood * Snedecor's F * Gamma distribution * more gamma stuff * Beroulli() * inverse Gaussian * Getting modifed GLM.jl tests to pass * drop pdf and logpdf * Poisson bugfix * update Normal(μ, σ²) * Gamma(μ, ϕ) for GLMs * updates for GLM support * start on truncated * update parameterized measures * drop FactoredBase * drop old LazyArrays dependency * insupport(::Distribution) * Left out"Dists." * don't export `ifelse` (https://github.com/cscherrer/MeasureTheory.jl/issues/192) * Kleisli => TransitionKernel * depend on StatsBase * tests passing * bump MeasureBase version * work on truncated and censored * improve func_string * Simplify logdensity_def(::For, x) * Move truncated and censored updates to separate branches * newline * comment out in-progress stuff * newline * bump version * update formatting spec * more formatting * tweedie docs * drop redundant exports * update exports * omega => lambda * drop SequentialEx * get tests passing * add kernel tests * gitignore * better `Pretty.tile` for Affine and AffineTransforms * formatting * kleisli => kernel * update tile(::For) * update Compat version * bump MB version * update gamma * Let's come back to InverseGaussian * CI on 1.7 * update IfElse * formatting --- .JuliaFormatter.toml | 2 +- .github/workflows/ci.yml | 1 + .gitignore | 8 ++ Project.toml | 15 ++-- docs/src/affine.md | 10 +-- src/MeasureTheory.jl | 32 ++++---- src/combinators/affine.jl | 76 ++++++++++--------- src/combinators/chain.jl | 3 +- src/combinators/exponential-families.jl | 31 +++++--- src/combinators/for.jl | 78 +++++-------------- src/combinators/product.jl | 14 ++-- src/combinators/transforms.jl | 22 ++++-- src/combinators/tweedie.jl | 99 +++++++++++++++++++------ src/combinators/weighted.jl | 2 +- src/const.jl | 2 +- src/distproxy.jl | 24 ++++-- src/distributions.jl | 4 +- src/macros.jl | 1 + src/parameterized.jl | 2 +- src/parameterized/bernoulli.jl | 26 +++++-- src/parameterized/beta.jl | 8 +- src/parameterized/binomial.jl | 22 +++++- src/parameterized/cauchy.jl | 9 ++- src/parameterized/dirichlet.jl | 2 +- src/parameterized/exponential.jl | 2 +- src/parameterized/gamma.jl | 67 +++++++++++++++++ src/parameterized/gumbel.jl | 2 +- src/parameterized/inverse-gamma.jl | 5 +- src/parameterized/inverse-gaussian.jl | 53 +++++++++++++ src/parameterized/laplace.jl | 9 ++- src/parameterized/lkj-cholesky.jl | 6 +- src/parameterized/multinomial.jl | 3 +- src/parameterized/mvnormal.jl | 12 +-- src/parameterized/negativebinomial.jl | 2 + src/parameterized/normal.jl | 25 ++++--- src/parameterized/poisson.jl | 18 +++-- src/parameterized/snedecorf.jl | 34 +++++++++ src/parameterized/studentt.jl | 13 ++-- src/realized.jl | 9 ++- src/smart-constructors.jl | 11 +-- src/transforms/ordered.jl | 2 +- src/utils.jl | 3 +- test/runtests.jl | 65 ++++++++++------ 43 files changed, 563 insertions(+), 271 deletions(-) create mode 100644 src/parameterized/gamma.jl create mode 100644 src/parameterized/inverse-gaussian.jl create mode 100644 src/parameterized/snedecorf.jl diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml index 154894d0..df2313d9 100644 --- a/.JuliaFormatter.toml +++ b/.JuliaFormatter.toml @@ -6,7 +6,7 @@ whitespace_ops_in_indices = false remove_extra_newlines = true import_to_using = false pipe_to_function_call = false -short_to_long_function_def = false +short_to_long_function_def = true always_use_return = false whitespace_in_kwargs = true annotate_untyped_fields_with_any = false diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ce867052..193fa517 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -23,6 +23,7 @@ jobs: matrix: version: - '1.6' + - '1.7' os: - ubuntu-latest - macOS-latest diff --git a/.gitignore b/.gitignore index a89f99e4..8a1bd58e 100644 --- a/.gitignore +++ b/.gitignore @@ -9,3 +9,11 @@ /.vscode/ /coverage coverage-lcov.info +paper/paper.aux +paper/paper.bbl +paper/paper.blg +paper/paper.fdb_latexmk +paper/paper.fls +paper/paper.log +paper/paper.out +paper/paper.synctex.gz diff --git a/Project.toml b/Project.toml index 5567b347..79caed99 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "MeasureTheory" uuid = "eadaa1a4-d27c-401d-8699-e962e1bbc33b" authors = ["Chad Scherrer and contributors"] -version = "0.15.1" +version = "0.16.0" [deps] Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697" @@ -12,14 +12,12 @@ ConstructionBase = "187b0558-2788-49d3-abe0-74a17ed4e7c9" DensityInterface = "b429d917-457f-4dbc-8f4c-0cc954292b1d" Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" DynamicIterators = "6c76993d-992e-5bf1-9e63-34920a5a5a38" -FLoops = "cc61a311-1640-44b5-9fba-1b764f453329" FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b" IfElse = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173" Infinities = "e1ba4f0e-776d-440f-acd9-e1d2e9742647" InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240" InverseFunctions = "3587e190-3f89-42d0-90ee-14403ec27112" KeywordCalls = "4d827475-d3e4-43d6-abe3-9688362ede9f" -LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" @@ -35,6 +33,8 @@ Reexport = "189a3867-3050-52da-a836-e630ba90ab69" SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b" Static = "aedffcd0-7271-4cad-89d0-dc628f76c6d3" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c" TransformVariables = "84d833dd-6860-57f9-a1a7-6da5db126cff" Tricks = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775" @@ -42,24 +42,22 @@ Tricks = "410a4b4d-49e4-4fbc-ab6d-cb71b17b3775" [compat] Accessors = "0.1" ChangesOfVariables = "0.1" -Compat = "3.35" +Compat = "3.35, 4" ConcreteStructs = "0.2" ConstructionBase = "1.3" DensityInterface = "0.4" Distributions = "0.25" DynamicIterators = "0.4" -FLoops = "0.2" FillArrays = "0.12, 0.13" IfElse = "0.1" Infinities = "0.1" InverseFunctions = "0.1" KeywordCalls = "0.2" -LazyArrays = "0.21, 0.22" LogExpFunctions = "0.3.3" MLStyle = "0.4" MacroTools = "0.5" MappedArrays = "0.4" -MeasureBase = "0.7" +MeasureBase = "0.9" NamedTupleTools = "0.13, 0.14" NestedTuples = "0.3" PositiveFactorizations = "0.2" @@ -68,7 +66,8 @@ Reexport = "1" SpecialFunctions = "1, 2" Static = "0.5, 0.6" StaticArrays = "1.3" -StatsFuns = "0.9" +StatsBase = "0.32, 0.33" +StatsFuns = "0.9, 1" TransformVariables = "0.5, 0.6" Tricks = "0.1" julia = "1.6" diff --git a/docs/src/affine.md b/docs/src/affine.md index 00975861..40293af8 100644 --- a/docs/src/affine.md +++ b/docs/src/affine.md @@ -69,17 +69,17 @@ Here the `- logdet(σ)` is the "log absolute Jacobian", required to account for The above requires solving a linear system, which adds some overhead. Even with the convenience of a lower triangular system, it's still not quite as efficient as multiplication. -In addition to the covariance ``Σ``, it's also common to parameterize a multivariate normal by its _precision matrix_, defined as the inverse of the covariance matrix, ``Ω = Σ⁻¹``. Similar to our use of ``σ`` for the lower Cholesky factor of `Σ`, we'll use ``ω`` for the lower Cholesky factor of ``Ω``. +In addition to the covariance ``Σ``, it's also common to parameterize a multivariate normal by its _precision matrix_, defined as the inverse of the covariance matrix, ``Ω = Σ⁻¹``. Similar to our use of ``σ`` for the lower Cholesky factor of `Σ`, we'll use ``λ`` for the lower Cholesky factor of ``Ω``. This parameterization enables more efficient calculation of the log-density using only multiplication and addition, ```julia -logdensity_def(d::Normal{(:μ,:ω)}, x) = logdensity_def(d.ω * (x - d.μ)) + logdet(d.ω) +logdensity_def(d::Normal{(:μ,:λ)}, x) = logdensity_def(d.λ * (x - d.μ)) + logdet(d.λ) ``` ## `AffineTransform` -Transforms like ``z → σ z + μ`` and ``z → ω \ z + μ`` can be specified in MeasureTheory.jl using an `AffineTransform`. For example, +Transforms like ``z → σ z + μ`` and ``z → λ \ z + μ`` can be specified in MeasureTheory.jl using an `AffineTransform`. For example, ```julia julia> f = AffineTransform((μ=3.,σ=2.)) @@ -91,11 +91,11 @@ julia> f(1.0) In the univariate case this is relatively simple to invert. But if `σ` is a matrix, matrix inversion becomes necessary. This is not always possible as lower triangular matrices are not closed under matrix inversion and as such are not guaranteed to exist. -With multiple parameterizations of a given family of measures, we can work around these issues. The inverse transform of a ``(μ,σ)`` transform will be in terms of ``(μ,ω)``, and vice-versa. So +With multiple parameterizations of a given family of measures, we can work around these issues. The inverse transform of a ``(μ,σ)`` transform will be in terms of ``(μ,λ)``, and vice-versa. So ```julia julia> f⁻¹ = inverse(f) -AffineTransform{(:μ, :ω), Tuple{Float64, Float64}}((μ = -1.5, ω = 2.0)) +AffineTransform{(:μ, :λ), Tuple{Float64, Float64}}((μ = -1.5, λ = 2.0)) julia> f(f⁻¹(4)) 4.0 diff --git a/src/MeasureTheory.jl b/src/MeasureTheory.jl index da97a78e..de66967a 100644 --- a/src/MeasureTheory.jl +++ b/src/MeasureTheory.jl @@ -1,7 +1,6 @@ module MeasureTheory using Random -using FLoops using MeasureBase using MLStyle @@ -9,13 +8,14 @@ using NestedTuples import TransformVariables const TV = TransformVariables -using TransformVariables: as, asℝ₊, as𝕀, asℝ +using TransformVariables: asℝ₊, as𝕀, asℝ, transform import Base import Distributions const Dists = Distributions export TV +export transform export ≪ export gentype export For @@ -30,7 +30,7 @@ export CountingMeasure export TrivialMeasure export Likelihood export testvalue -export basekleisli +export basekernel using Infinities using DynamicIterators @@ -45,7 +45,6 @@ import LogExpFunctions import NamedTupleTools import InverseFunctions: inverse export inverse -export ifelse import MeasureBase: insupport, instance, marginals import MeasureBase: @@ -53,7 +52,7 @@ import MeasureBase: logdensity_def, density_def, basemeasure, - kleisli, + kernel, params, paramnames, ∫, @@ -62,7 +61,13 @@ import MeasureBase: import MeasureBase: ≪ using MeasureBase: BoundedInts, BoundedReals, CountingMeasure, IntegerDomain, IntegerNumbers using MeasureBase: weightedmeasure, restrict -using MeasureBase: AbstractKleisli +using MeasureBase: AbstractTransitionKernel + +import Statistics: mean, var, std + +import MeasureBase: likelihood +export likelihood +export log_likelihood_ratio using StaticArrays @@ -74,8 +79,8 @@ import Base: rand using Reexport @reexport using MeasureBase -import IfElse: ifelse -@reexport using IfElse +import IfElse +using IfElse using Tricks: static_hasmethod @@ -87,6 +92,7 @@ export AffineTransform export insupport export For +using MeasureBase: kernel using MeasureBase: Returns import MeasureBase: proxy, @useproxy import MeasureBase: basemeasure_depth @@ -101,10 +107,6 @@ gentype(μ::AbstractMeasure) = typeof(testvalue(μ)) # gentype(μ::AbstractMeasure) = gentype(basemeasure(μ)) -import Distributions: logpdf, pdf - -export pdf, logpdf - xlogx(x::Number) = LogExpFunctions.xlogx(x) xlogx(x, y) = x * log(x) @@ -114,6 +116,8 @@ xlogy(x, y) = x * log(y) xlog1py(x::Number, y::Number) = LogExpFunctions.xlog1py(x, y) xlog1py(x, y) = x * log(1 + y) +as(args...; kwargs...) = TV.as(args...; kwargs...) + include("utils.jl") include("const.jl") include("combinators/for.jl") @@ -126,7 +130,6 @@ include("combinators/weighted.jl") include("combinators/product.jl") include("combinators/transforms.jl") include("combinators/exponential-families.jl") - include("resettable-rng.jl") include("realized.jl") include("combinators/chain.jl") @@ -151,6 +154,9 @@ include("parameterized/binomial.jl") include("parameterized/multinomial.jl") include("parameterized/lkj-cholesky.jl") include("parameterized/negativebinomial.jl") +include("parameterized/gamma.jl") +include("parameterized/snedecorf.jl") +include("parameterized/inverse-gaussian.jl") include("combinators/ifelse.jl") diff --git a/src/combinators/affine.jl b/src/combinators/affine.jl index 298dcd33..c705bd02 100644 --- a/src/combinators/affine.jl +++ b/src/combinators/affine.jl @@ -3,9 +3,9 @@ using LinearAlgebra const AFFINEPARS = [ (:μ, :σ) - (:μ, :ω) + (:μ, :λ) (:σ,) - (:ω,) + (:λ,) (:μ,) ] @@ -13,7 +13,13 @@ struct AffineTransform{N,T} par::NamedTuple{N,T} end -quoteof(f::AffineTransform) = :(AffineTransform($(quoteof(f.par)))) +function Pretty.tile(f::AffineTransform) + result = Pretty.literal("AffineTransform") + result *= Pretty.literal(sprint(show, params(f); context = :compact => true)) + result +end + +Base.show(io::IO, f::AffineTransform) = Pretty.pprint(io, f) params(f::AffineTransform) = getfield(f, :par) @@ -23,17 +29,17 @@ Base.propertynames(d::AffineTransform{N}) where {N} = N import InverseFunctions: inverse -@inline inverse(f::AffineTransform{(:μ, :σ)}) = AffineTransform((μ = -(f.σ \ f.μ), ω = f.σ)) -@inline inverse(f::AffineTransform{(:μ, :ω)}) = AffineTransform((μ = -f.ω * f.μ, σ = f.ω)) -@inline inverse(f::AffineTransform{(:σ,)}) = AffineTransform((ω = f.σ,)) -@inline inverse(f::AffineTransform{(:ω,)}) = AffineTransform((σ = f.ω,)) +@inline inverse(f::AffineTransform{(:μ, :σ)}) = AffineTransform((μ = -(f.σ \ f.μ), λ = f.σ)) +@inline inverse(f::AffineTransform{(:μ, :λ)}) = AffineTransform((μ = -f.λ * f.μ, σ = f.λ)) +@inline inverse(f::AffineTransform{(:σ,)}) = AffineTransform((λ = f.σ,)) +@inline inverse(f::AffineTransform{(:λ,)}) = AffineTransform((σ = f.λ,)) @inline inverse(f::AffineTransform{(:μ,)}) = AffineTransform((μ = -f.μ,)) # `size(f) == (m,n)` means `f : ℝⁿ → ℝᵐ` Base.size(f::AffineTransform{(:μ, :σ)}) = size(f.σ) -Base.size(f::AffineTransform{(:μ, :ω)}) = size(f.ω) +Base.size(f::AffineTransform{(:μ, :λ)}) = size(f.λ) Base.size(f::AffineTransform{(:σ,)}) = size(f.σ) -Base.size(f::AffineTransform{(:ω,)}) = size(f.ω) +Base.size(f::AffineTransform{(:λ,)}) = size(f.λ) function Base.size(f::AffineTransform{(:μ,)}) (n,) = size(f.μ) @@ -44,9 +50,9 @@ Base.size(f::AffineTransform, n::Int) = @inbounds size(f)[n] (f::AffineTransform{(:μ,)})(x) = x + f.μ (f::AffineTransform{(:σ,)})(x) = f.σ * x -(f::AffineTransform{(:ω,)})(x) = f.ω \ x +(f::AffineTransform{(:λ,)})(x) = f.λ \ x (f::AffineTransform{(:μ, :σ)})(x) = f.σ * x + f.μ -(f::AffineTransform{(:μ, :ω)})(x) = f.ω \ x + f.μ +(f::AffineTransform{(:μ, :λ)})(x) = f.λ \ x + f.μ rowsize(x) = () rowsize(x::AbstractArray) = (size(x, 1),) @@ -78,13 +84,13 @@ end return x end -@inline function apply!(x, f::AffineTransform{(:ω,),Tuple{F}}, z) where {F<:Factorization} - ldiv!(x, f.ω, z) +@inline function apply!(x, f::AffineTransform{(:λ,),Tuple{F}}, z) where {F<:Factorization} + ldiv!(x, f.λ, z) return x end -@inline function apply!(x, f::AffineTransform{(:ω,)}, z) - ldiv!(x, factorize(f.ω), z) +@inline function apply!(x, f::AffineTransform{(:λ,)}, z) + ldiv!(x, factorize(f.λ), z) return x end @@ -94,8 +100,8 @@ end return x end -@inline function apply!(x, f::AffineTransform{(:μ, :ω)}, z) - apply!(x, AffineTransform((ω = f.ω,)), z) +@inline function apply!(x, f::AffineTransform{(:μ, :λ)}, z) + apply!(x, AffineTransform((λ = f.λ,)), z) apply!(x, AffineTransform((μ = f.μ,)), x) return x end @@ -113,9 +119,9 @@ logjac(x::Number) = log(abs(x)) # TODO: `log` doesn't work for the multivariate case, we need the log absolute determinant logjac(f::AffineTransform{(:μ, :σ)}) = logjac(f.σ) -logjac(f::AffineTransform{(:μ, :ω)}) = -logjac(f.ω) +logjac(f::AffineTransform{(:μ, :λ)}) = -logjac(f.λ) logjac(f::AffineTransform{(:σ,)}) = logjac(f.σ) -logjac(f::AffineTransform{(:ω,)}) = -logjac(f.ω) +logjac(f::AffineTransform{(:λ,)}) = -logjac(f.λ) logjac(f::AffineTransform{(:μ,)}) = 0.0 ############################################################################### @@ -136,7 +142,9 @@ struct Affine{N,M,T} <: AbstractMeasure end function Pretty.tile(d::Affine) - Pretty.list_layout(Pretty.tile.([params(d.f), d.parent]); prefix = :Affine) + pars = Pretty.literal(sprint(show, params(d.f); context = :compact => true)) + + Pretty.list_layout([pars, Pretty.tile(d.parent)]; prefix = :Affine) end function testvalue(d::Affine) @@ -175,15 +183,15 @@ end Base.size(d::Affine) = size(d.μ) Base.size(d::Affine{(:σ,)}) = (size(d.σ, 1),) -Base.size(d::Affine{(:ω,)}) = (size(d.ω, 2),) +Base.size(d::Affine{(:λ,)}) = (size(d.λ, 2),) @inline function logdensity_def(d::Affine{(:σ,)}, x::AbstractArray) z = solve(d.σ, x) MeasureBase.unsafe_logdensityof(d.parent, z) end -@inline function logdensity_def(d::Affine{(:ω,)}, x::AbstractArray) - z = d.ω * x +@inline function logdensity_def(d::Affine{(:λ,)}, x::AbstractArray) + z = d.λ * x MeasureBase.unsafe_logdensityof(d.parent, z) end @@ -197,8 +205,8 @@ end MeasureBase.unsafe_logdensityof(d.parent, z) end -@inline function logdensity_def(d::Affine{(:μ, :ω)}, x::AbstractArray) - z = d.ω * mappedarray(-, x, d.μ) +@inline function logdensity_def(d::Affine{(:μ, :λ)}, x::AbstractArray) + z = d.λ * mappedarray(-, x, d.μ) MeasureBase.unsafe_logdensityof(d.parent, z) end @@ -207,7 +215,7 @@ end logdensity_def(d.parent, z) end -# # # logdensity_def(d::Affine{(:μ,:ω)}, x) = logdensity_def(d.parent, d.σ \ (x - d.μ)) +# # # logdensity_def(d::Affine{(:μ,:λ)}, x) = logdensity_def(d.parent, d.σ \ (x - d.μ)) # # @inline function logdensity_def(d::Affine{(:μ,:σ), P, Tuple{V,M}}, x) where {P, V<:AbstractVector, M<:AbstractMatrix} # # z = x - d.μ # # σ = d.σ @@ -219,10 +227,10 @@ end # # sum(zⱼ -> logdensity_def(d.parent, zⱼ), z) # # end -# # # logdensity_def(d::Affine{(:μ,:ω)}, x) = logdensity_def(d.parent, d.ω * (x - d.μ)) -# # @inline function logdensity_def(d::Affine{(:μ,:ω), P,Tuple{V,M}}, x) where {P,V<:AbstractVector, M<:AbstractMatrix} +# # # logdensity_def(d::Affine{(:μ,:λ)}, x) = logdensity_def(d.parent, d.λ * (x - d.μ)) +# # @inline function logdensity_def(d::Affine{(:μ,:λ), P,Tuple{V,M}}, x) where {P,V<:AbstractVector, M<:AbstractMatrix} # # z = x - d.μ -# # lmul!(d.ω, z) +# # lmul!(d.λ, z) # # logdensity_def(d.parent, z) # # end @@ -240,10 +248,12 @@ end # We can't do this until we know we're working with Lebesgue measure, since for # example it wouldn't make sense to apply a log-Jacobian to a point measure -@inline basemeasure(d::Affine{N,L}) where {N,L<:Lebesgue} = +@inline function basemeasure(d::Affine{N,L}) where {N,L<:Lebesgue} weightedmeasure(-logjac(d), d.parent) -@inline basemeasure(d::Affine{N,L}) where {N,L<:LebesgueMeasure} = +end +@inline function basemeasure(d::Affine{N,L}) where {N,L<:LebesgueMeasure} weightedmeasure(-logjac(d), d.parent) +end logjac(d::Affine) = logjac(getfield(d, :f)) @@ -267,9 +277,9 @@ end # end supportdim(nt::NamedTuple{(:μ, :σ)}) = colsize(nt.σ) -supportdim(nt::NamedTuple{(:μ, :ω)}) = rowsize(nt.ω) +supportdim(nt::NamedTuple{(:μ, :λ)}) = rowsize(nt.λ) supportdim(nt::NamedTuple{(:σ,)}) = colsize(nt.σ) -supportdim(nt::NamedTuple{(:ω,)}) = rowsize(nt.ω) +supportdim(nt::NamedTuple{(:λ,)}) = rowsize(nt.λ) supportdim(nt::NamedTuple{(:μ,)}) = size(nt.μ) asparams(::Affine, ::StaticSymbol{:μ}) = asℝ diff --git a/src/combinators/chain.jl b/src/combinators/chain.jl index de5803c5..187c3227 100644 --- a/src/combinators/chain.jl +++ b/src/combinators/chain.jl @@ -63,8 +63,9 @@ struct DynamicFor{T,K,S} <: AbstractMeasure iter::S end -Pretty.quoteof(r::DynamicFor) = +function Pretty.quoteof(r::DynamicFor) :(DynamicFor($(Pretty.quoteof(r.κ)), $(Pretty.quoteof(r.iter)))) +end function DynamicFor(κ::K, iter::S) where {K,S} T = typeof(κ(first(iter))) diff --git a/src/combinators/exponential-families.jl b/src/combinators/exponential-families.jl index 68bd930e..525fde4f 100644 --- a/src/combinators/exponential-families.jl +++ b/src/combinators/exponential-families.jl @@ -1,6 +1,7 @@ export ExponentialFamily -@concrete terse struct ExponentialFamily <: AbstractKleisli +@concrete terse struct ExponentialFamily <: AbstractTransitionKernel + support_contains base mdim pdim @@ -9,15 +10,26 @@ export ExponentialFamily a end -function ExponentialFamily(base, mdim, pdim, t, a) - return ExponentialFamily(base, mdim, pdim, t, I, a) +MeasureBase.insupport(fam::ExponentialFamily, x) = fam.support_contains(x) + +function ExponentialFamily(support_contains, base, mdim, pdim, t, a) + return ExponentialFamily(support_contains, base, mdim, pdim, t, I, a) end function MeasureBase.powermeasure(fam::ExponentialFamily, dims::NTuple{N,I}) where {N,I} + support_contains(x) = all(xj -> fam.support_contains(xj), x) t = Tuple((y -> f.(y) for f in fam.t)) a(η) = BroadcastArray(fam.a, η) p = prod(dims) - ExponentialFamily(fam.base^dims, fam.mdim * p, fam.pdim * p, t, fam.x, a) + ExponentialFamily( + support_contains, + fam.base^dims, + fam.mdim * p, + fam.pdim * p, + t, + fam.x, + a, + ) end @concrete terse struct ExpFamMeasure <: AbstractMeasure @@ -26,6 +38,8 @@ end a # instantiated to a value end +MeasureBase.insupport(μ::ExpFamMeasure, x) = μ.fam.support_contains(x) + @inline function (fam::ExponentialFamily)(β) η = fam.x * β a = fam.a(η) @@ -76,19 +90,18 @@ end c end -export likelihood +# function regression(fam, uᵀ, vᵀ) -function regression(fam, uᵀ, vᵀ) -end +# end -function likelihood(fam::ExponentialFamily, y) +function MeasureBase.likelihood(fam::ExponentialFamily, y) c = logdensityof(fam.base, y) t = ApplyArray(vcat, (f.(y) for f in fam.t)...) tᵀx = t' * fam.x ExpFamLikelihood(fam, y, tᵀx, c) end -@inline function logdensity_def(ℓ::ExpFamLikelihood, β) +@inline function logdensityof(ℓ::ExpFamLikelihood, β) xβ = ApplyArray(*, ℓ.fam.x, β) a = sum(ℓ.fam.a(xβ)) # a = sum(ℓ.fam.a, ApplyArray(*, ℓ.fam.uᵀ', ℓ.fam.vᵀ, β)) diff --git a/src/combinators/for.jl b/src/combinators/for.jl index dec16fa4..1227cbf6 100644 --- a/src/combinators/for.jl +++ b/src/combinators/for.jl @@ -12,16 +12,6 @@ struct For{T,F,I} <: AbstractProductMeasure end @inline For{T,F,I}(f::F, inds::I) where {T,F,I} = new{T,F,I}(f, inds) - - function For{Union{},F,I}(f::F, inds::I) where {F,I} - println("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX") - @warn "Empty `For` construction. This should not be happening" - @show f(first(zip(inds...))...) - println.(stacktrace()) - println("XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX") - # @error "Empty `For` construction" - return ProductMeasure(mappedarray(i -> f(Tuple(i)...), CartesianIndices(inds...))) - end end @generated function For(f::F, inds::I) where {F,I<:Tuple} @@ -33,6 +23,8 @@ end end end +as(d::For) = as(Array, as(first(marginals(d))), size(first(d.inds))...) + # For(f, gen::Base.Generator) = ProductMeasure(Base.Generator(f ∘ gen.f, gen.iter)) # @inline function tailcall_iter(f, iter) @@ -69,54 +61,22 @@ end @inline function logdensity_def( d::For{T,F,I}, x::AbstractVector{X}; - exec = SequentialEx(simd = true), ) where {X,T,F,I<:Tuple{<:AbstractVector}} ind = only(d.inds) - @floop exec for j in eachindex(x) - local i = getindex(ind, j) - local Δℓ = @inbounds logdensity_def(d.f(i), x[j]) - @reduce ℓ += Δℓ + sum(eachindex(x)) do j + i = getindex(ind, j) + @inbounds logdensity_def(d.f(i), x[j]) end - ℓ end -function logdensity_def(d::For, x::AbstractVector; exec = SequentialEx(simd = true)) - @floop exec for j in eachindex(x) - local i = (getindex(ind, j) for ind in d.inds) - local Δℓ = @inbounds logdensity_def(d.f(i...), x[j]) - @reduce ℓ += Δℓ +function logdensity_def(d::For, x::AbstractVector) + get_i(j) = tuple((getindex(ind, j) for ind in d.inds)...) + sum(eachindex(x)) do j + i = get_i(j) + @inbounds logdensity_def(d.f(i...), x[j]) end - - ℓ -end - -function logdensity_def( - d::For{T,F,I}, - x::AbstractArray{X}; - exec = SequentialEx(simd = true), -) where {T,F,I,X} - @floop exec for j in CartesianIndices(x) - local i = (getindex(ind, j) for ind in d.inds) - local Δℓ = @inbounds logdensity_def(d.f(i...), x[j]) - @reduce ℓ += Δℓ - end - - ℓ end -# function logdensity_def(d::For{T,F,I}, x) where {N,T,F,I<:NTuple{N,<:Base.Generator}} -# sum(zip(x, d.inds...)) do (xⱼ, dⱼ...) -# logdensity_def(d.f(dⱼ...), xⱼ) -# end -# end - -# function logdensity_def(d::For{T,F,I}, x::AbstractVector) where {N,T,F,I<:NTuple{N,<:Base.Generator}} - -# sum(zip(x, d.inds...)) do (xⱼ, dⱼ...) -# logdensity_def(d.f(dⱼ...), xⱼ) -# end -# end - function marginals(d::For{T,F,Tuple{I}}) where {T,F,I} f = d.f data = first(d.inds) @@ -136,13 +96,10 @@ function marginals(d::For{T,F,I}) where {N,T,F,I<:NTuple{N,<:Base.Generator}} Iterators.map(d.f, d.inds...) end -@generated function basemeasure(d::For{T,F,I}) where {T,F,I} - B = Core.Compiler.return_type(basemeasure, Tuple{T}) +function basemeasure(d::For{T,F,I}) where {T,F,I} + B = typeof(basemeasure(d.f(map(first, d.inds)...))) sing = static(Base.issingletontype(B)) - quote - $(Expr(:meta, :inline)) - _basemeasure(d, $B, $sing) - end + _basemeasure(d, B, sing) end @inline function _basemeasure( @@ -175,7 +132,7 @@ end ::Type{B}, ::False, ) where {T,F,I,B<:AbstractMeasure} - f = basekleisli(d.f) + f = basekernel(d.f) For{B}(f, d.inds) end @@ -196,7 +153,7 @@ function _basemeasure( ::Type{B}, ::False, ) where {N,T<:AbstractMeasure,F,I<:NTuple{N,<:Base.Generator},B} - f = basekleisli(d.f) + f = basekernel(d.f) For{B}(f, d.inds) end @@ -206,7 +163,7 @@ function Pretty.tile(d::For{T}) where {T} result *= Pretty.literal("}") result *= Pretty.list_layout([ Pretty.literal(func_string(d.f, Tuple{_eltype.(d.inds)...})), - Pretty.tile.(d.inds)..., + Pretty.literal(sprint(show, d.inds; context = :compact => true)), ]) end @@ -283,8 +240,9 @@ julia> For(eachrow(rand(4,2))) do x Normal(x[1], x[2]) end |> marginals |> colle """ @inline For{T}(f, inds...) where {T} = For{T}(f, inds) @inline For{T}(f, n::Integer) where {T} = For{T}(f, static(1):n) -@inline For{T}(f, inds::Integer...) where {T} = +@inline function For{T}(f, inds::Integer...) where {T} For{T}(i -> f(Tuple(i)...), Base.CartesianIndices(inds)) +end @inline For(f, inds...) = For(f, inds) @inline For(f, n::Integer) = For(f, static(1):n) diff --git a/src/combinators/product.jl b/src/combinators/product.jl index 19dd2112..dc5ffd9c 100644 --- a/src/combinators/product.jl +++ b/src/combinators/product.jl @@ -1,21 +1,25 @@ -function TV.as(d::PowerMeasure) +function as(d::PowerMeasure) as(Array, as(d.parent), length.(d.axes)...) end -function TV.as(d::ProductMeasure{A}) where {A<:AbstractArray} +function as(d::ProductMeasure{A}) where {A<:AbstractArray} d1 = marginals(d).f(first(marginals(d).data)) - as(Array, TV.as(d1), size(marginals(d))...) + as(Array, as(d1), size(marginals(d))...) +end + +function as(d::ProductMeasure{T}) where {T<:Tuple} + as(map(as, d.marginals)) end ############################################################################### # I <: Base.Generator -function TV.as(d::ProductMeasure{<:Base.Generator}) +function as(d::ProductMeasure{<:Base.Generator}) d1 = marginals(d).f(first(marginals(d).iter)) as(Array, as(d1), size(marginals(d))...) end -# function TV.as(d::ProductMeasure{Returns{T},F,A}) where {T,F,A<:AbstractArray} +# function as(d::ProductMeasure{Returns{T},F,A}) where {T,F,A<:AbstractArray} # as(Array, as(d.f.f.value), size(d.xs)) # end diff --git a/src/combinators/transforms.jl b/src/combinators/transforms.jl index 11279803..24a39957 100644 --- a/src/combinators/transforms.jl +++ b/src/combinators/transforms.jl @@ -72,8 +72,9 @@ end Pullback(f::AbstractTransform, ν, logjac = True()) = Pullback(TV.transform(f), ν, logjac) -Pushforward(f::AbstractTransform, ν, logjac = True()) = +function Pushforward(f::AbstractTransform, ν, logjac = True()) Pushforward(TV.transform(f), ν, logjac) +end Pullback(f::CallableInverse, ν, logjac = True()) = Pushforward(TV.transform(f.t), ν, logjac) @@ -91,24 +92,29 @@ basemeasure(μ::Pullback) = Pullback(μ.f, basemeasure(μ.ν), False()) basemeasure(ν::Pushforward) = Pushforward(ν.f, basemeasure(ν.μ), False()) -TV.as(ν::Pushforward) = ν.f ∘ as(ν.μ) +as(ν::Pushforward) = ν.f ∘ as(ν.μ) -TV.as(μ::Pullback) = TV.inverse(μ.f) ∘ μ.ν +as(μ::Pullback) = TV.inverse(μ.f) ∘ μ.ν -TV.as(::Lebesgue) = asℝ +as(::Lebesgue) = asℝ # TODO: Make this work for affine embeddings -TV.as(d::Affine) = _as_affine(_firstval(d)) +as(d::Affine) = _as_affine(_firstval(d)) _firstval(d::Affine) = first(values(getfield(getfield(d, :f), :par))) _as_affine(x::Real) = asℝ _as_affine(x::AbstractArray) = as(Vector, size(x, 1)) -basemeasure( +function basemeasure( ::Pushforward{TV.CallableTransform{T},Lebesgue{ℝ}}, -) where {T<:TV.ScalarTransform} = Lebesgue(ℝ) -basemeasure(::Pullback{TV.CallableTransform{T},Lebesgue{ℝ}}) where {T<:TV.ScalarTransform} = +) where {T<:TV.ScalarTransform} + Lebesgue(ℝ) +end +function basemeasure( + ::Pullback{TV.CallableTransform{T},Lebesgue{ℝ}}, +) where {T<:TV.ScalarTransform} Lebesgue(ℝ) +end # t = as𝕀 # μ = Normal() # ν = Pushforward(t, μ) diff --git a/src/combinators/tweedie.jl b/src/combinators/tweedie.jl index 14ad2ed1..12b4318f 100644 --- a/src/combinators/tweedie.jl +++ b/src/combinators/tweedie.jl @@ -1,37 +1,59 @@ export Tweedie -struct Tweedie{B,Θ,P} <: AbstractKleisli +abstract type AbstractEDM <: AbstractTransitionKernel end + +""" +From https://en.wikipedia.org/wiki/Tweedie_distribution: + +> The Tweedie distributions include a number of familiar distributions as well as +some unusual ones, each being specified by the domain of the index parameter. We +have the + + extreme stable distribution, p < 0, + normal distribution, p = 0, + Poisson distribution, p = 1, + compound Poisson-gamma distribution, 1 < p < 2, + gamma distribution, p = 2, + positive stable distributions, 2 < p < 3, + Inverse Gaussian distribution, p = 3, + positive stable distributions, p > 3, and + extreme stable distributions, p = ∞. + +For 0 < p < 1 no Tweedie model exists. Note that all stable distributions mean +actually generated by stable distributions. +""" +struct Tweedie{S,B,D,P} <: AbstractEDM + support_contains::S base::B - θ::Θ + dim::D p::P end struct TweedieMeasure{B,Θ,P,S,C} <: AbstractMeasure - base::B + fam::Tweedie{B,Θ,P} θ::Θ - p::P σ::S cumulant::C end -@inline function tweedie_mean(::StaticFloat64{0.0}, θ) - return θ -end +mean(d::TweedieMeasure) = tweedie_mean(d.fam.p, d.θ) -@inline function tweedie_mean(::StaticFloat64{1.0}, θ) - return exp(θ) -end - -@inline function tweedie_mean(::StaticFloat64{2.0}, θ) - return inv(log(-θ)) -end +var(d::TweedieMeasure) = d.σ^2 * mean(d)^d.fam.p -@generated function tweedie_mean(::StaticFloat64{p}, θ) where {p} - α_minus_1 = (p - 2) / (p - 1) - 1 +############################################################################### +# Tweedie cumulants - quote - $(Expr(:meta, :inline)) - (θ / α_minus_1)^α_minus_1 +@inline function tweedie_cumulant(p::P, θ) where {P} + if p == zero(P) + return 0.5 * θ^2 + elseif p == one(P) + return exp(θ) + elseif p == 2 + return -log(-θ) + else + α = (p - 2) / (p - 1) + coeff = (α - 1) / α + return coeff * (θ / (α - 1))^α end end @@ -66,6 +88,43 @@ end TweedieMeasure(base, θ, p, σ) end +############################################################################### +# Tweedie mean function + +@inline function tweedie_mean(p::P, θ) where {P} + if p == zero(P) + return θ + elseif p == one(P) + return exp(θ) + elseif p == 2 + return inv(log(-θ)) + else + α_minus_1 = (p - 2) / (p - 1) - 1 + return (θ / α_minus_1)^α_minus_1 + end +end + +@inline function tweedie_mean(::StaticFloat64{0.0}, θ) + return θ +end + +@inline function tweedie_mean(::StaticFloat64{1.0}, θ) + return exp(θ) +end + +@inline function tweedie_mean(::StaticFloat64{2.0}, θ) + return inv(log(-θ)) +end + +@generated function tweedie_mean(::StaticFloat64{p}, θ) where {p} + α_minus_1 = (p - 2) / (p - 1) - 1 + + quote + $(Expr(:meta, :inline)) + (θ / α_minus_1)^α_minus_1 + end +end + basemeasure(d::TweedieMeasure) = d.base function logdensity_def(d::TweedieMeasure, x) @@ -86,8 +145,6 @@ struct TweedieLikelihood{C,Θ,H,T,A} <: AbstractLikelihood a::A end -export likelihood - function likelihood(fam::Tweedie, x) c = logdensityof(fam.base, x) t = fam.t(x) diff --git a/src/combinators/weighted.jl b/src/combinators/weighted.jl index 84aae8c4..e7357644 100644 --- a/src/combinators/weighted.jl +++ b/src/combinators/weighted.jl @@ -1,2 +1,2 @@ -TV.as(μ::AbstractWeightedMeasure) = TV.as(μ.base) +as(μ::AbstractWeightedMeasure) = as(μ.base) diff --git a/src/const.jl b/src/const.jl index 5f4e2add..3bfd6fb0 100644 --- a/src/const.jl +++ b/src/const.jl @@ -5,7 +5,7 @@ end asConst(x) = AsConst(x) -TV.as(d::Dirac) = AsConst(d.x) +as(d::Dirac) = AsConst(d.x) TV.dimension(t::AsConst) = 0 diff --git a/src/distproxy.jl b/src/distproxy.jl index d1e17c0f..60701217 100644 --- a/src/distproxy.jl +++ b/src/distproxy.jl @@ -1,12 +1,24 @@ export proxy function proxy end -PROXIES = Dict(:Distributions => [ - :mean - :std - :entropy - :cdf -]) +import Statistics +import StatsBase + +PROXIES = Dict( + :Statistics => [ + :mean + :std + :var + :quantile + ], + :StatsBase => [:entropy], + :Distributions => [ + :cdf + :ccdf + :logcdf + :logccdf + ], +) for m in keys(PROXIES) for f in PROXIES[m] diff --git a/src/distributions.jl b/src/distributions.jl index a8aab4d2..09fa90ee 100644 --- a/src/distributions.jl +++ b/src/distributions.jl @@ -38,6 +38,8 @@ end logdensity_def(μ::Dists.Distribution, x) = Dists.logpdf(μ, x) -density(μ::Dists.Distribution, x) = Dists.pdf(μ, x) +density_def(μ::Dists.Distribution, x) = Dists.pdf(μ, x) testvalue(d::Dists.Distribution) = rand(d) + +insupport(d::Dists.Distribution, x) = logdensityof(d, x) > -Inf diff --git a/src/macros.jl b/src/macros.jl index 998645aa..50e5b7e0 100644 --- a/src/macros.jl +++ b/src/macros.jl @@ -148,6 +148,7 @@ function _half(__module__, ex) halfdist = esc(Symbol(:Half, dist)) dist = esc(dist) quote + $halfdist(arg::NamedTuple) = Half($dist(arg)) $halfdist(args...) = Half($dist(args...)) $halfdist(; kwargs...) = Half($dist(; kwargs...)) end diff --git a/src/parameterized.jl b/src/parameterized.jl index 354f7d01..751989df 100644 --- a/src/parameterized.jl +++ b/src/parameterized.jl @@ -74,4 +74,4 @@ function asparams(μ::M, nt::NamedTuple = NamedTuple()) where {M<:ParameterizedM asparams(constructorof(M), nt) end -TV.as(::Half) = asℝ₊ +as(::Half) = asℝ₊ diff --git a/src/parameterized/bernoulli.jl b/src/parameterized/bernoulli.jl index 7f83a25c..10db3c46 100644 --- a/src/parameterized/bernoulli.jl +++ b/src/parameterized/bernoulli.jl @@ -3,10 +3,14 @@ export Bernoulli import Base -@parameterized Bernoulli(p) +@parameterized Bernoulli() + +@kwstruct Bernoulli() @kwstruct Bernoulli(p) +Bernoulli(p) = Bernoulli((p = p,)) + basemeasure(::Bernoulli) = CountingMeasure() testvalue(::Bernoulli) = true @@ -15,11 +19,18 @@ insupport(::Bernoulli, x) = x == true || x == false @inline function logdensity_def(d::Bernoulli{(:p,)}, y) p = d.p - f = ifelse(y, () -> log(p), () -> log(1 - p)) - return f() + y == true ? log(p) : log1p(-p) +end + +function density_def(::Bernoulli{()}, y) + return 0.5 end -function density(d::Bernoulli{(:p,)}, y) +@inline function logdensity_def(d::Bernoulli{()}, y) + return -logtwo +end + +function density_def(d::Bernoulli{(:p,)}, y) p = d.p return 2 * p * y - p - y + 1 end @@ -29,17 +40,20 @@ end return y * x - log1pexp(x) end -function density(d::Bernoulli{(:logitp,)}, y) +function density_def(d::Bernoulli{(:logitp,)}, y) exp_x = exp(d.logitp) return exp_x^y / (1 + exp_x) end gentype(::Bernoulli) = Bool +Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{()}) = rand(rng, T) < one(T) / 2 + Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:p,)}) = rand(rng, T) < d.p -Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:logitp,)}) = +function Base.rand(rng::AbstractRNG, T::Type, d::Bernoulli{(:logitp,)}) rand(rng, T) < logistic(d.logitp) +end asparams(::Type{<:Bernoulli}, ::StaticSymbol{:p}) = as𝕀 asparams(::Type{<:Bernoulli}, ::StaticSymbol{:logitp}) = asℝ diff --git a/src/parameterized/beta.jl b/src/parameterized/beta.jl index 5a9d9dc9..1763261d 100644 --- a/src/parameterized/beta.jl +++ b/src/parameterized/beta.jl @@ -13,17 +13,15 @@ export Beta beta => β ] -TV.as(::Beta) = as𝕀 +as(::Beta) = as𝕀 @inline function logdensity_def(d::Beta{(:α, :β),Tuple{A,B}}, x::X) where {A,B,X} return xlogy(d.α - 1, x) + xlog1py(d.β - 1, -x) end @inline function basemeasure(d::Beta{(:α, :β)}) - constℓ = 0.0 - varℓ = Returns(-logbeta(d.α, d.β)) - base = Lebesgue(ℝ) - FactoredBase(constℓ, varℓ, base) + ℓ = -logbeta(d.α, d.β) + weightedmeasure(ℓ, Lebesgue()) end Base.rand(rng::AbstractRNG, T::Type, μ::Beta) = rand(rng, Dists.Beta(μ.α, μ.β)) diff --git a/src/parameterized/binomial.jl b/src/parameterized/binomial.jl index 44314be7..72981b96 100644 --- a/src/parameterized/binomial.jl +++ b/src/parameterized/binomial.jl @@ -13,6 +13,22 @@ basemeasure(d::Binomial) = CountingMeasure() testvalue(::Binomial) = 0 +@kwstruct Binomial() + +@inline function logdensity_def(d::Binomial{()}, y) + return -logtwo +end + +@inline function insupport(d::Binomial{()}, x) + x ∈ (0, 1) +end + +function Base.rand(rng::AbstractRNG, ::Type, d::Binomial{(:n, :p)}) + rand(rng, Dists.Binomial(d.n, d.p)) +end + +Binomial(n) = Binomial(n, 0.5) + ############################################################################### @kwstruct Binomial(n, p) @@ -68,10 +84,12 @@ function Base.rand( end proxy(d::Binomial{(:n, :p),Tuple{I,A}}) where {I<:Integer,A} = Dists.Binomial(d.n, d.p) -proxy(d::Binomial{(:n, :logitp),Tuple{I,A}}) where {I<:Integer,A} = +function proxy(d::Binomial{(:n, :logitp),Tuple{I,A}}) where {I<:Integer,A} Dists.Binomial(d.n, logistic(d.logitp)) -proxy(d::Binomial{(:n, :probitp),Tuple{I,A}}) where {I<:Integer,A} = +end +function proxy(d::Binomial{(:n, :probitp),Tuple{I,A}}) where {I<:Integer,A} Dists.Binomial(d.n, Φ(d.probitp)) +end asparams(::Type{<:Binomial}, ::StaticSymbol{:p}) = as𝕀 asparams(::Type{<:Binomial}, ::StaticSymbol{:logitp}) = asℝ diff --git a/src/parameterized/cauchy.jl b/src/parameterized/cauchy.jl index c5edf635..44331a4a 100644 --- a/src/parameterized/cauchy.jl +++ b/src/parameterized/cauchy.jl @@ -9,8 +9,8 @@ export Cauchy, HalfCauchy @kwstruct Cauchy(μ) @kwstruct Cauchy(σ) @kwstruct Cauchy(μ, σ) -@kwstruct Cauchy(ω) -@kwstruct Cauchy(μ, ω) +@kwstruct Cauchy(λ) +@kwstruct Cauchy(μ, λ) logdensity_def(d::Cauchy, x) = logdensity_def(proxy(d), x) @@ -32,14 +32,15 @@ for N in AFFINEPARS @eval begin proxy(d::Cauchy{$N}) = affine(params(d), Cauchy()) @useproxy Cauchy{$N} - rand(rng::AbstractRNG, ::Type{T}, d::Cauchy{$N}) where {T} = + function rand(rng::AbstractRNG, ::Type{T}, d::Cauchy{$N}) where {T} rand(rng, T, affine(params(d), Cauchy())) + end end end ≪(::Cauchy, ::Lebesgue{X}) where {X<:Real} = true -TV.as(::Cauchy) = asℝ +as(::Cauchy) = asℝ @half Cauchy diff --git a/src/parameterized/dirichlet.jl b/src/parameterized/dirichlet.jl index 853164b5..bd2940b0 100644 --- a/src/parameterized/dirichlet.jl +++ b/src/parameterized/dirichlet.jl @@ -7,7 +7,7 @@ using FillArrays @parameterized Dirichlet(α) -TV.as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α)) +as(d::Dirichlet{(:α,)}) = TV.UnitSimplex(length(d.α)) @inline function basemeasure(μ::Dirichlet{(:α,)}) α = μ.α diff --git a/src/parameterized/exponential.jl b/src/parameterized/exponential.jl index e7cb2347..09e5a2e9 100644 --- a/src/parameterized/exponential.jl +++ b/src/parameterized/exponential.jl @@ -16,7 +16,7 @@ end Base.rand(rng::AbstractRNG, T::Type, μ::Exponential{()}) = randexp(rng, T) -TV.as(::Exponential) = asℝ₊ +as(::Exponential) = asℝ₊ ########################## # Scale β diff --git a/src/parameterized/gamma.jl b/src/parameterized/gamma.jl new file mode 100644 index 00000000..73681efc --- /dev/null +++ b/src/parameterized/gamma.jl @@ -0,0 +1,67 @@ + +# Gamma distribution + +export Gamma + +@parameterized Gamma() + +@kwstruct Gamma() + +proxy(::Gamma{()}) = Exponential() + +@useproxy Gamma{()} + +@kwstruct Gamma(k) + +@inline function logdensity_def(d::Gamma{(:k,)}, x) + return xlogy(d.k - 1, x) - x +end + +function basemeasure(d::Gamma{(:k,)}) + ℓ = -loggamma(d.k) + weightedmeasure(ℓ, Lebesgue()) +end + +@kwstruct Gamma(k, σ) + +Gamma(k, σ) = Gamma((k = k, σ = σ)) + +@useproxy Gamma{(:k, :σ)} +function proxy(d::Gamma{(:k, :σ)}) + affine(NamedTuple{(:σ,)}(d.σ), Gamma((k = d.k,))) +end + +@kwstruct Gamma(k, λ) + +@useproxy Gamma{(:k, :λ)} +function proxy(d::Gamma{(:k, :λ)}) + affine(NamedTuple{(:λ,)}(d.λ), Gamma((k = d.k,))) +end + +Base.rand(rng::AbstractRNG, T::Type, μ::Gamma{()}) = rand(rng, T, Exponential()) + +Base.rand(rng::AbstractRNG, T::Type, μ::Gamma{(:k,)}) = rand(rng, Dists.Gamma(μ.k)) + +as(::Gamma) = asℝ₊ + +insupport(::Gamma, x) = x > 0 + +@kwstruct Gamma(μ, ϕ) + +@inline function logdensity_def(d::Gamma{(:μ, :ϕ)}, x) + x_μ = x / d.μ + inv(d.ϕ) * (log(x_μ) - x_μ) - log(x) +end + +function basemeasure(d::Gamma{(:μ, :ϕ)}) + ϕ = d.ϕ + ϕinv = inv(ϕ) + ℓ = -ϕinv * log(ϕ) - first(logabsgamma(ϕinv)) + weightedmeasure(ℓ, Lebesgue()) +end + +function basemeasure(d::Gamma{(:μ, :ϕ),Tuple{M,StaticFloat64{ϕ}}}) where {M,ϕ} + ϕinv = inv(ϕ) + ℓ = static(-ϕinv * log(ϕ) - first(logabsgamma(ϕinv))) + weightedmeasure(ℓ, Lebesgue()) +end diff --git a/src/parameterized/gumbel.jl b/src/parameterized/gumbel.jl index 2bf1cd30..a3a6394c 100644 --- a/src/parameterized/gumbel.jl +++ b/src/parameterized/gumbel.jl @@ -33,7 +33,7 @@ function Base.rand(rng::AbstractRNG, d::Gumbel{()}) -log(-log(u)) end -TV.as(::Gumbel) = asℝ +as(::Gumbel) = asℝ ≪(::Gumbel, ::Lebesgue{X}) where {X<:Real} = true diff --git a/src/parameterized/inverse-gamma.jl b/src/parameterized/inverse-gamma.jl index f2e5b87c..c488e35e 100644 --- a/src/parameterized/inverse-gamma.jl +++ b/src/parameterized/inverse-gamma.jl @@ -12,11 +12,12 @@ export InverseGamma return xlogy(α + 1, xinv) - xinv - loggamma(α) end -Base.rand(rng::AbstractRNG, T::Type, μ::InverseGamma{(:shape,)}) = +function Base.rand(rng::AbstractRNG, T::Type, μ::InverseGamma{(:shape,)}) rand(rng, Dists.InverseGamma(μ.shape)) +end ≪(::InverseGamma, ::Lebesgue{X}) where {X<:Real} = true -TV.as(::InverseGamma) = asℝ₊ +as(::InverseGamma) = asℝ₊ # @μσ_methods InverseGamma(shape) diff --git a/src/parameterized/inverse-gaussian.jl b/src/parameterized/inverse-gaussian.jl new file mode 100644 index 00000000..382fecbf --- /dev/null +++ b/src/parameterized/inverse-gaussian.jl @@ -0,0 +1,53 @@ +# InverseGaussian distribution + +export InverseGaussian + +@parameterized InverseGaussian() + +# @kwstruct InverseGaussian() +# @kwstruct InverseGaussian(μ) +# @kwstruct InverseGaussian(μ, λ) + +# InverseGaussian(μ) = InverseGaussian((μ = μ,)) + +# InverseGaussian(μ, λ) = InverseGaussian((μ = μ, λ = λ)) + +# @useproxy InverseGaussian{(:k, :θ)} + +# function logdensity_def(d::InverseGaussian{()}, x) +# return (-3log(x) - (x - 1)^2 / x) / 2 +# end + +# function logdensity_def(d::InverseGaussian{(:μ,)}, x) +# μ = d.μ +# return (-3log(x) - (x - μ)^2 / (μ^2 * x)) / 2 +# end + +# function logdensity_def(d::InverseGaussian{(:μ, :λ)}, x) +# μ, λ = d.μ, d.λ +# return (log(λ) - 3log(x) - λ * (x - μ)^2 / (μ^2 * x)) / 2 +# end + +# function basemeasure(::InverseGaussian) +# ℓ = static(-log2π / 2) +# weightedmeasure(ℓ, Lebesgue()) +# end + +Base.rand(rng::AbstractRNG, T::Type, d::InverseGaussian) = rand(rng, proxy(d)) + +as(::InverseGaussian) = asℝ₊ + +insupport(::InverseGaussian, x) = x > 0 + +# GLM parameterization +@kwstruct InverseGaussian(μ, ϕ) + +function logdensity_def(d::InverseGaussian{(:μ, :ϕ)}, x) + μ, ϕ = d.μ, d.ϕ + return (-3log(x) - (x - μ)^2 / (ϕ * μ^2 * x)) / 2 +end + +function basemeasure(d::InverseGaussian{(:μ, :ϕ)}) + ℓ = static(-0.5) * (static(log2π) + log(d.ϕ)) + weightedmeasure(ℓ, Lebesgue()) +end diff --git a/src/parameterized/laplace.jl b/src/parameterized/laplace.jl index 2471f326..ae913bb4 100644 --- a/src/parameterized/laplace.jl +++ b/src/parameterized/laplace.jl @@ -9,8 +9,8 @@ export Laplace @kwstruct Laplace(μ) @kwstruct Laplace(σ) @kwstruct Laplace(μ, σ) -@kwstruct Laplace(ω) -@kwstruct Laplace(μ, ω) +@kwstruct Laplace(λ) +@kwstruct Laplace(μ, λ) for N in AFFINEPARS @eval begin @@ -31,10 +31,11 @@ basemeasure(::Laplace{()}) = WeightedMeasure(static(-logtwo), Lebesgue(ℝ)) # @affinepars Laplace -Base.rand(rng::AbstractRNG, ::Type{T}, μ::Laplace{()}) where {T} = +function Base.rand(rng::AbstractRNG, ::Type{T}, μ::Laplace{()}) where {T} rand(rng, Dists.Laplace()) +end Base.rand(rng::AbstractRNG, ::Type{T}, μ::Laplace) where {T} = Base.rand(rng, T, proxy(μ)) ≪(::Laplace, ::Lebesgue{X}) where {X<:Real} = true -TV.as(::Laplace) = asℝ +as(::Laplace) = asℝ diff --git a/src/parameterized/lkj-cholesky.jl b/src/parameterized/lkj-cholesky.jl index 65bd7296..1133f50c 100644 --- a/src/parameterized/lkj-cholesky.jl +++ b/src/parameterized/lkj-cholesky.jl @@ -77,16 +77,16 @@ end return s end -TV.as(d::LKJCholesky) = CorrCholesky(d.k) +as(d::LKJCholesky) = CorrCholesky(d.k) @inline function basemeasure(d::LKJCholesky{(:k, :η)}) - t = TV.as(d) + t = as(d) base = Pushforward(t, Lebesgue(ℝ)^TV.dimension(t), False()) WeightedMeasure(Dists.lkj_logc0(d.k, d.η), base) end @inline function basemeasure(d::LKJCholesky{(:k, :logη)}) - t = TV.as(d) + t = as(d) η = exp(d.logη) base = Pushforward(t, Lebesgue(ℝ)^TV.dimension(t), False()) WeightedMeasure(Dists.lkj_logc0(d.k, η), base) diff --git a/src/parameterized/multinomial.jl b/src/parameterized/multinomial.jl index 6ececb92..781de7d7 100644 --- a/src/parameterized/multinomial.jl +++ b/src/parameterized/multinomial.jl @@ -24,8 +24,9 @@ end return s end -Base.rand(rng::AbstractRNG, T::Type, μ::Multinomial) = +function Base.rand(rng::AbstractRNG, T::Type, μ::Multinomial) rand(rng, Dists.Multinomial(μ.n, μ.p)) +end proxy(d::Multinomial{(:p,)}) = Dists.Multinomial(d.n, d.p) diff --git a/src/parameterized/mvnormal.jl b/src/parameterized/mvnormal.jl index 51e3a414..21c38ec5 100644 --- a/src/parameterized/mvnormal.jl +++ b/src/parameterized/mvnormal.jl @@ -9,9 +9,9 @@ export MvNormal @kwstruct MvNormal(μ) @kwstruct MvNormal(σ) -@kwstruct MvNormal(ω) +@kwstruct MvNormal(λ) @kwstruct MvNormal(μ, σ) -@kwstruct MvNormal(μ, ω) +@kwstruct MvNormal(μ, λ) supportdim(d::MvNormal) = supportdim(params(d)) @@ -33,8 +33,8 @@ insupport(::MvNormal, x) = true # affine(nt, Normal() ^ dim) # end -# function MvNormal(nt::NamedTuple{(:ω,)}) -# dim = rowsize(nt.ω) +# function MvNormal(nt::NamedTuple{(:λ,)}) +# dim = rowsize(nt.λ) # affine(nt, Normal() ^ dim) # end @@ -43,7 +43,7 @@ insupport(::MvNormal, x) = true # affine(nt, Normal() ^ dim) # end -# function MvNormal(nt::NamedTuple{(:μ, :ω,)}) -# dim = rowsize(nt.ω) +# function MvNormal(nt::NamedTuple{(:μ, :λ,)}) +# dim = rowsize(nt.λ) # affine(nt, Normal() ^ dim) # end diff --git a/src/parameterized/negativebinomial.jl b/src/parameterized/negativebinomial.jl index cf860e24..dcca6741 100644 --- a/src/parameterized/negativebinomial.jl +++ b/src/parameterized/negativebinomial.jl @@ -18,6 +18,8 @@ end ############################################################################### @kwstruct NegativeBinomial(r, p) +NegativeBinomial(n) = NegativeBinomial(n, 0.5) + @inline function logdensity_def(d::NegativeBinomial{(:r, :p)}, y) (r, p) = (d.r, d.p) return -log(y + r) - logbeta(r, y + 1) + xlogy(y, p) + xlog1py(r, -p) diff --git a/src/parameterized/normal.jl b/src/parameterized/normal.jl index cd0cfd62..03e2e537 100644 --- a/src/parameterized/normal.jl +++ b/src/parameterized/normal.jl @@ -42,8 +42,8 @@ insupport(d::Normal) = Returns(true) @kwstruct Normal(μ) @kwstruct Normal(σ) @kwstruct Normal(μ, σ) -@kwstruct Normal(ω) -@kwstruct Normal(μ, ω) +@kwstruct Normal(λ) +@kwstruct Normal(μ, λ) params(::Type{N}) where {N<:Normal} = () @@ -51,7 +51,7 @@ Normal(μ, σ) = Normal((μ = μ, σ = σ)) Normal(nt::NamedTuple{N,Tuple{Vararg{AbstractArray}}}) where {N} = MvNormal(nt) -TV.as(::Normal) = asℝ +as(::Normal) = asℝ # `@kwalias` defines some alias names, giving users flexibility in the names # they use. For example, σ² is standard notation for the variance parameter, but @@ -63,6 +63,7 @@ TV.as(::Normal) = asℝ std => σ sigma => σ var => σ² + ϕ => σ² ] # It's often useful to be able to map into the parameter space for a given @@ -134,22 +135,24 @@ Base.rand(rng::Random.AbstractRNG, ::Type{T}, μ::Normal) where {T} = rand(rng, @half Normal # A single unnamed parameter for `HalfNormal` should be interpreted as a `σ` -HalfNormal(σ) = HalfNormal(σ = σ) +HalfNormal(σ) = HalfNormal((σ = σ,)) ############################################################################### +@kwstruct Normal(σ²) @kwstruct Normal(μ, σ²) -@inline function logdensity_def(d::Normal{(:σ²)}, x) - σ² = d.σ² - -0.5 * (log(σ²) + (x^2 / σ²)) +@inline function logdensity_def(d::Normal{(:σ²,)}, x) + -0.5 * x^2 / d.σ² end -@inline function logdensity_def(d::Normal{(:μ, :σ²)}, x) - μ = d.μ - σ² = d.σ² - -0.5 * (log(σ²) + ((x - μ)^2 / σ²)) +@inline function basemeasure(d::Normal{(:σ²,)}) + ℓ = static(-0.5) * (static(log2π) + log(d.σ²)) + weightedmeasure(ℓ, Lebesgue()) end +proxy(d::Normal{(:μ, :σ²)}) = affine((μ = d.μ,), Normal((σ² = d.σ²,))) +@useproxy Normal{(:μ, :σ²)} + ############################################################################### @kwstruct Normal(μ, τ) diff --git a/src/parameterized/poisson.jl b/src/parameterized/poisson.jl index 9191692c..7135fd70 100644 --- a/src/parameterized/poisson.jl +++ b/src/parameterized/poisson.jl @@ -2,23 +2,30 @@ export Poisson import Base -using SpecialFunctions: logfactorial +using SpecialFunctions: loggamma -@parameterized Poisson(λ) +@parameterized Poisson() +@kwstruct Poisson() @kwstruct Poisson(λ) +Poisson(λ) = Poisson((λ = λ,)) + basemeasure(::Poisson) = Counting(BoundedInts(static(0), static(Inf))) +@inline function logdensity_def(d::Poisson{()}, y::T) where {T} + return y - one(T) - loggamma(one(T) + y) +end + @inline function logdensity_def(d::Poisson{(:λ,)}, y) λ = d.λ - return y * log(λ) - λ - logfactorial(y) + return y * log(λ) - λ - loggamma(1 + y) end @kwstruct Poisson(logλ) @inline function logdensity_def(d::Poisson{(:logλ,)}, y) - return y * d.logλ - exp(d.logλ) - logfactorial(y) + return y * d.logλ - exp(d.logλ) - loggamma(1 + y) end asparams(::Type{<:Poisson}, ::StaticSymbol{:λ}) = asℝ₊ @@ -27,8 +34,9 @@ asparams(::Type{<:Poisson}, ::StaticSymbol{:logλ}) = asℝ gentype(::Poisson) = Int Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:λ,)}) = rand(rng, Dists.Poisson(d.λ)) -Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:logλ,)}) = +function Base.rand(rng::AbstractRNG, T::Type, d::Poisson{(:logλ,)}) rand(rng, Dists.Poisson(exp(d.logλ))) +end @inline function insupport(::Poisson, x) isinteger(x) && x ≥ 0 diff --git a/src/parameterized/snedecorf.jl b/src/parameterized/snedecorf.jl new file mode 100644 index 00000000..cdce60fb --- /dev/null +++ b/src/parameterized/snedecorf.jl @@ -0,0 +1,34 @@ + +# Snedecor's F distribution + +export SnedecorF + +@parameterized SnedecorF(ν1, ν2) + +@kwstruct SnedecorF(ν1, ν2) + +@inline function logdensity_def(d::SnedecorF{(:ν1, :ν2)}, y) + ν1, ν2 = d.ν1, d.ν2 + ν1ν2 = ν1 / ν2 + val = (xlogy(ν1, ν1ν2) + xlogy(ν1 - 2, y) - xlogy(ν1 + ν2, 1 + ν1ν2 * y)) / 2 + return val +end + +@inline function basemeasure(d::SnedecorF{(:ν1, :ν2)}) + ℓ = -logbeta(d.ν1 / 2, d.ν2 / 2) + weightedmeasure(ℓ, Lebesgue()) +end + +xform(::SnedecorF) = asℝ₊ + +Base.rand(rng::AbstractRNG, T::Type, d::SnedecorF) = rand(rng, proxy(d)) + +proxy(d::SnedecorF{(:ν1, :ν2)}) = Dists.FDist(d.ν1, d.ν2) + +asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν1}) = asℝ₊ +asparams(::Type{<:SnedecorF}, ::StaticSymbol{:ν2}) = asℝ₊ + +insupport(::SnedecorF, x) = x > 0 + +# cdf(d::SnedecorF, x) = StatsFuns.fdistcdf(d.ν1, d.ν2, x) +# ccdf(d::SnedecorF, x) = StatsFuns.fdistccdf(d.ν1, d.ν2, x) diff --git a/src/parameterized/studentt.jl b/src/parameterized/studentt.jl index 12ded3fb..cf08fef7 100644 --- a/src/parameterized/studentt.jl +++ b/src/parameterized/studentt.jl @@ -8,14 +8,15 @@ export StudentT, HalfStudentT @kwstruct StudentT(ν) @kwstruct StudentT(ν, μ) @kwstruct StudentT(ν, σ) -@kwstruct StudentT(ν, ω) +@kwstruct StudentT(ν, λ) @kwstruct StudentT(ν, μ, σ) -@kwstruct StudentT(ν, μ, ω) +@kwstruct StudentT(ν, μ, λ) for N in AFFINEPARS @eval begin - proxy(d::StudentT{(:ν, $N...)}) = + function proxy(d::StudentT{(:ν, $N...)}) affine(NamedTuple{$N}(params(d)), StudentT((ν = d.ν,))) + end end end @@ -48,10 +49,8 @@ end end @inline function basemeasure(d::StudentT{(:ν,)}) - constℓ = 0.0 - varℓ() = loggamma((d.ν + 1) / 2) - loggamma(d.ν / 2) - log(π * d.ν) / 2 - base = Lebesgue(ℝ) - FactoredBase(constℓ, varℓ, base) + ℓ = loggamma((d.ν + 1) / 2) - loggamma(d.ν / 2) - log(π * d.ν) / 2 + weightedmeasure(ℓ, Lebesgue(ℝ)) end xform(::StudentT) = asℝ diff --git a/src/realized.jl b/src/realized.jl index 0eb9b413..cda405d8 100644 --- a/src/realized.jl +++ b/src/realized.jl @@ -32,8 +32,9 @@ reset!(rv::Realized) = reset!(rv.rng) Base.show(io::IO, r::Realized) = Pretty.pprint(io, r) -Pretty.quoteof(r::Realized) = +function Pretty.quoteof(r::Realized) :(Realized($(Pretty.quoteof(r.rng)), $(Pretty.quoteof(r.iter)))) +end Base.IteratorEltype(mc::Realized) = Base.HasEltype() @@ -96,8 +97,9 @@ end Base.show(io::IO, r::RealizedMeasures) = Pretty.pprint(io, r) -Pretty.quoteof(r::RealizedMeasures) = +function Pretty.quoteof(r::RealizedMeasures) :(RealizedMeasures($(Pretty.quoteof(r.rng)), $(Pretty.quoteof(r.iter)))) +end function Base.iterate(rm::RealizedMeasures, s = nothing) val, s = iterate(rm, s) @@ -119,8 +121,9 @@ RealizedSamples(rng::AbstractRNG, iter) = RealizedSamples(Realized(rng, iter)) Base.show(io::IO, r::RealizedSamples) = Pretty.pprint(io, r) -Pretty.quoteof(r::RealizedSamples) = +function Pretty.quoteof(r::RealizedSamples) :(RealizedSamples($(Pretty.quoteof(r.parent.rng)), $(Pretty.quoteof(r.parent.iter)))) +end function iter_helper(x) isnothing(x) && return x diff --git a/src/smart-constructors.jl b/src/smart-constructors.jl index 7de14054..97219470 100644 --- a/src/smart-constructors.jl +++ b/src/smart-constructors.jl @@ -11,19 +11,12 @@ function affine(f::AffineTransform, parent::WeightedMeasure) WeightedMeasure(parent.logweight, affine(f, parent.base)) end -function affine(f::AffineTransform, parent::FactoredBase) - constℓ = parent.constℓ - varℓ = parent.varℓ - base = affine(f, parent.base) - FactoredBase(constℓ, varℓ, base) -end - using MeasureBase: RealNumbers function affine(f::AffineTransform{(:μ, :σ)}, parent::Lebesgue{RealNumbers}) affine(AffineTransform((σ = f.σ,)), parent) end -function affine(f::AffineTransform{(:μ, :ω)}, parent::Lebesgue{RealNumbers}) - affine(AffineTransform((ω = f.ω,)), parent) +function affine(f::AffineTransform{(:μ, :λ)}, parent::Lebesgue{RealNumbers}) + affine(AffineTransform((λ = f.λ,)), parent) end diff --git a/src/transforms/ordered.jl b/src/transforms/ordered.jl index 6b63fe07..7c927a66 100644 --- a/src/transforms/ordered.jl +++ b/src/transforms/ordered.jl @@ -37,7 +37,7 @@ function TV.transform_with(flag::TV.LogJacFlag, t::Ordered, x, index::T) where { x = mappedarray(xj -> xj + OrderedΔx, x) - @inbounds (y[1], ℓ, _) = TV.transform_with(flag, TV.as(Real, lo, hi), x, index) + @inbounds (y[1], ℓ, _) = TV.transform_with(flag, as(Real, lo, hi), x, index) index += 1 @inbounds for i in 2:len diff --git a/src/utils.jl b/src/utils.jl index cca327c7..ba2e2b7c 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -1,6 +1,5 @@ using LinearAlgebra export mydot -using LazyArrays function solve(A::Union{AbstractMatrix,Factorization}, y::AbstractArray) (m, n) = size(A) @@ -89,7 +88,7 @@ function func_string(f, types) result *= string(s[1]) result *= "->" else - result *= string(tuple(s...)) + result *= "(" * join(string.(tuple(s...)), ", ") * ")" result *= "->" end c1 = string(c[1]) diff --git a/test/runtests.jl b/test/runtests.jl index 3fe5f7a2..8959532d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -9,6 +9,7 @@ using FillArrays using MeasureTheory using MeasureBase.Interface +using MeasureTheory: kernel using Aqua Aqua.test_all(MeasureTheory; ambiguities = false, unbound_args = false) @@ -179,9 +180,25 @@ end end end -@testset "Kleisli" begin - κ = kleisli(Dirac) +@testset "TransitionKernel" begin + κ = kernel(Dirac) @test rand(κ(1.1)) == 1.1 + + k1 = kernel() do x + Normal(x, x^2) + end + + k2 = kernel(Normal) do x + (μ = x, σ = x^2) + end + + k3 = kernel(Normal; μ = identity, σ = abs2) + + k4 = kernel(Normal; μ = first, σ = last) do x + (x, x^2) + end + + @test k1(3) == k2(3) == k3(3) == k4(3) == Normal(3, 9) end @testset "For" begin @@ -209,10 +226,10 @@ end end end -import MeasureTheory.:⋅ +import MeasureTheory: ⋅ -function ⋅(μ::Normal, kleisli) - m = kleisli(μ) +function ⋅(μ::Normal, kernel) + m = kernel(μ) Normal(μ = m.μ.μ, σ = sqrt(m.μ.σ^2 + m.σ^2)) end struct AffineMap{S,T} @@ -250,11 +267,11 @@ end # Q = 0.2 # μ = Normal(μ=ξ0, σ=sqrt(P0)) -# kleisli = MeasureTheory.kleisli(Normal; μ=AffineMap(Φ, β), σ=MeasureTheory.AsConst(Q)) +# kernel = MeasureTheory.kernel(Normal; μ=AffineMap(Φ, β), σ=MeasureTheory.AsConst(Q)) -# @test (μ ⋅ kleisli).μ == Normal(μ = 0.9, σ = 0.824621).μ +# @test (μ ⋅ kernel).μ == Normal(μ = 0.9, σ = 0.824621).μ -# chain = Chain(kleisli, μ) +# chain = Chain(kernel, μ) # dyniterate(iter::TimeLift, ::Nothing) = dyniterate(iter, 0=>nothing) # tr1 = trace(TimeLift(chain), nothing, u -> u[1] > 15) @@ -292,12 +309,12 @@ end ℓs = [ Likelihood(Normal{(:μ,)}, 3.0), - # Likelihood(kleisli(Normal, x -> (μ=x, σ=2.0)), 3.0) + # Likelihood(kernel(Normal, x -> (μ=x, σ=2.0)), 3.0) ] for (d, p) in dps for ℓ in ℓs - @test logdensity_def(d ⊙ ℓ, p) ≈ logdensity_def(d, p) + logdensity_def(ℓ, p) + @test logdensityof(d ⊙ ℓ, p) ≈ logdensityof(d, p) + logdensityof(ℓ.k(p), ℓ.x) end end end @@ -499,7 +516,7 @@ end end for (M, nt) in testmeasures - for p in [(μ = 1,), (μ = 1, σ = 2), (μ = 1, ω = 2), (σ = 2,), (ω = 2,)] + for p in [(μ = 1,), (μ = 1, σ = 2), (μ = 1, λ = 2), (σ = 2,), (λ = 2,)] d = M(merge(nt, p)) @info "Testing $d" test_noerrors(d) @@ -508,7 +525,7 @@ end # @show n # for k in 1:n # @show k - # pars = [(μ=randn(n),), (μ=randn(n),σ=randn(n,k)), (μ=randn(n),ω=randn(k,n)), (σ=randn(n,k),), (ω=randn(k,n),)] + # pars = [(μ=randn(n),), (μ=randn(n),σ=randn(n,k)), (μ=randn(n),λ=randn(k,n)), (σ=randn(n,k),), (λ=randn(k,n),)] # for p in pars # @show p # d = M(merge(nt, p)) @@ -524,7 +541,7 @@ end @test f(inverse(f)(1)) == 1 @test inverse(f)(f(1)) == 1 - f = AffineTransform((μ = 3, ω = 2)) + f = AffineTransform((μ = 3, λ = 2)) @test f(inverse(f)(1)) == 1 @test inverse(f)(f(1)) == 1 @@ -532,7 +549,7 @@ end @test f(inverse(f)(1)) == 1 @test inverse(f)(f(1)) == 1 - f = AffineTransform((ω = 2,)) + f = AffineTransform((λ = 2,)) @test f(inverse(f)(1)) == 1 @test inverse(f)(f(1)) == 1 @@ -544,10 +561,10 @@ end @testset "Affine" begin unif = ∫(x -> 0 < x < 1, Lebesgue(ℝ)) f1 = AffineTransform((μ = 3.0, σ = 2.0)) - f2 = AffineTransform((μ = 3.0, ω = 2.0)) + f2 = AffineTransform((μ = 3.0, λ = 2.0)) f3 = AffineTransform((μ = 3.0,)) f4 = AffineTransform((σ = 2.0,)) - f5 = AffineTransform((ω = 2.0,)) + f5 = AffineTransform((λ = 2.0,)) for f in [f1, f2, f3, f4, f5] par = getfield(f, :par) @@ -562,14 +579,14 @@ end σ = let x = randn(10, 3) cholesky(x' * x).L end - ω = inv(σ) + λ = inv(σ) x = randn(3) @test logdensity_def(Affine((μ = μ, σ = σ), d^3), x) ≈ - logdensity_def(Affine((μ = μ, ω = ω), d^3), x) + logdensity_def(Affine((μ = μ, λ = λ), d^3), x) @test logdensity_def(Affine((σ = σ,), d^3), x) ≈ - logdensity_def(Affine((ω = ω,), d^3), x) + logdensity_def(Affine((λ = λ,), d^3), x) @test logdensity_def(Affine((μ = μ,), d^3), x) ≈ logdensity_def(d^3, x - μ) d = ∫exp(x -> -x^2, Lebesgue(ℝ)) @@ -579,7 +596,7 @@ end @test logdensityof(a, x) ≈ logdensityof(d, inverse(a.f)(x)[1]) @test logdensityof(a, a.f(y)) ≈ logdensityof(d^1, y) - b = Affine((ω = [1 0]'',), d^1) + b = Affine((λ = [1 0]'',), d^1) @test logdensityof(b, x) ≈ logdensityof(d, inverse(b.f)(x)[1]) @test logdensityof(b, b.f(y)) ≈ logdensityof(d^1, y) end @@ -587,8 +604,10 @@ end @testset "IfElseMeasure" begin p = rand() x = randn() - @test logdensityof(MeasureTheory.ifelse(Bernoulli(p), Normal(), Normal()), x) ≈ + @test logdensityof(MeasureTheory.IfElse.ifelse(Bernoulli(p), Normal(), Normal()), x) ≈ logdensityof(Normal(), x) - @test logdensityof(MeasureTheory.ifelse(Bernoulli(p), Normal(2, 3), Normal()), x) ≈ - logdensityof(p * Normal(2, 3) + (1 - p) * Normal(), x) + @test logdensityof( + MeasureTheory.IfElse.ifelse(Bernoulli(p), Normal(2, 3), Normal()), + x, + ) ≈ logdensityof(p * Normal(2, 3) + (1 - p) * Normal(), x) end