diff --git a/src/MeasureBase.jl b/src/MeasureBase.jl index e29c4ae9..eeea2ad8 100644 --- a/src/MeasureBase.jl +++ b/src/MeasureBase.jl @@ -37,9 +37,7 @@ using Static using Static: StaticInteger using FunctionChains -export ≪ export gentype -export rebase export AbstractMeasure @@ -135,10 +133,8 @@ include("combinators/product.jl") include("combinators/power.jl") include("combinators/spikemixture.jl") include("combinators/likelihood.jl") -include("combinators/pointwise.jl") include("combinators/restricted.jl") include("combinators/smart-constructors.jl") -include("combinators/powerweighted.jl") include("combinators/conditional.jl") include("standard/stdmeasure.jl") @@ -156,6 +152,8 @@ include("density-core.jl") include("interface.jl") +include("measure_operators.jl") + using .Interface end # module MeasureBase diff --git a/src/absolutecontinuity.jl b/src/absolutecontinuity.jl index 8062198c..c65aeaf3 100644 --- a/src/absolutecontinuity.jl +++ b/src/absolutecontinuity.jl @@ -54,3 +54,6 @@ # representative(μ) ≪ representative(ν) && return true # return false # end + +# ≪(::M, ::WeightedMeasure{R,M}) where {R,M} = true +# ≪(::WeightedMeasure{R,M}, ::M) where {R,M} = true diff --git a/src/combinators/bind.jl b/src/combinators/bind.jl index cc2022f2..465f2bf7 100644 --- a/src/combinators/bind.jl +++ b/src/combinators/bind.jl @@ -1,36 +1,45 @@ +""" + struct MeasureBase.Bind{M,K} <: AbstractMeasure + +Represents a monatic bind. User code should create instances of `Bind` +directly, but should call `mbind(k, μ)` instead. +""" struct Bind{M,K} <: AbstractMeasure - μ::M k::K + μ::M end -export ↣ +getdof(d::Bind) = NoDOF{typeof(d)}() + +function Base.rand(rng::AbstractRNG, ::Type{T}, d::Bind) where {T} + x = rand(rng, T, d.μ) + y = rand(rng, T, d.k(x)) + return y +end """ -If -- μ is an `AbstractMeasure` or satisfies the Measure interface, and -- k is a function taking values from the support of μ and returning a measure + mbind(k, μ)::AbstractMeasure + +Given -Then `μ ↣ k` is a measure, called a *monadic bind*. In a -probabilistic programming language like Soss.jl, this could be expressed as +- a measure μ +- a kernel function k that takes values from the support of μ and returns a + measure -Note that bind is usually written `>>=`, but this symbol is unavailable in Julia. +The *monadic bind* operation `mbind(k, μ)` returns is a new measure. +If `ν == mbind(k, μ)` and all measures involved are sampleable, then +samples from `rand(ν)` follow the same distribution as those from `rand(k(rand(μ)))`. + + +A monadic bind ofen written as `>>=` (e.g. in Haskell), but this symbol is +unavailable in Julia. ``` -bind = @model μ,k begin - x ~ μ - y ~ k(x) - return y +μ = StdExponential() +ν = mbind(μ) do scale + pushfwd(Base.Fix1(*, scale), StdNormal()) end ``` - -See also `bind` and `Bind` """ -↣(μ, k) = bind(μ, k) - -bind(μ, k) = Bind(μ, k) - -function Base.rand(rng::AbstractRNG, ::Type{T}, d::Bind) where {T} - x = rand(rng, T, d.μ) - y = rand(rng, T, d.k(x)) - return y -end +mbind(k, μ) = Bind(k, μ) +export mbind diff --git a/src/combinators/likelihood.jl b/src/combinators/likelihood.jl index 6dfd164f..b244fd0f 100644 --- a/src/combinators/likelihood.jl +++ b/src/combinators/likelihood.jl @@ -11,9 +11,9 @@ abstract type AbstractLikelihood end # insupport(ℓ::AbstractLikelihood, p) = insupport(ℓ.k(p), ℓ.x) @doc raw""" - Likelihood(k::AbstractTransitionKernel, x) + Likelihood(k, x) -"Observe" a value `x`, yielding a function from the parameters to ℝ. +Default result of [`likelihoodof(k, x)`](@ref). Likelihoods are most commonly used in conjunction with an existing _prior_ measure to yield a new measure, the _posterior_. In Bayes's Law, we have @@ -64,39 +64,12 @@ With several parameters, things work as expected: --------- - Likelihood(M<:ParameterizedMeasure, constraint::NamedTuple, x) - -In some cases the measure might have several parameters, and we may want the -(log-)likelihood with respect to some subset of them. In this case, we can use -the three-argument form, where the second argument is a constraint. For example, - - julia> ℓ = Likelihood(Normal{(:μ,:σ)}, (σ=3.0,), 2.0) - Likelihood(Normal{(:μ, :σ), T} where T, (σ = 3.0,), 2.0) - -Similarly to the above, we have - - julia> density_def(ℓ, (μ=2.0,)) - 0.3333333333333333 - - julia> logdensity_def(ℓ, (μ=2.0,)) - -1.0986122886681098 - - julia> density_def(ℓ, 2.0) - 0.3333333333333333 - - julia> logdensity_def(ℓ, 2.0) - -1.0986122886681098 - ------------------------ - Finally, let's return to the expression for Bayes's Law, -``P(θ|x) ∝ P(θ) P(x|θ)`` +``P(θ|x) ∝ P(x|θ) P(θ)`` -The product on the right side is computed pointwise. To work with this in -MeasureBase, we have a "pointwise product" `⊙`, which takes a measure and a -likelihood, and returns a new measure, that is, the unnormalized posterior that -has density ``P(θ) P(x|θ)`` with respect to the base measure of the prior. +In measure theory, the product on the right side is the Lebesgue integral +of the likelihood with respect to the prior. For example, say we have @@ -104,23 +77,27 @@ For example, say we have x ~ Normal(μ,σ) σ = 1 -and we observe `x=3`. We can compute the posterior measure on `μ` as +and we observe `x=3`. We can compute the (non-normalized) posterior measure on +`μ` as - julia> post = Normal() ⊙ Likelihood(Normal{(:μ, :σ)}, (σ=1,), 3) - Normal() ⊙ Likelihood(Normal{(:μ, :σ), T} where T, (σ = 1,), 3) - - julia> logdensity_def(post, 2) - -2.5 + julia> prior = Normal() + julia> likelihood = Likelihood(μ -> Normal(μ, 1), 3) + julia> post = mintegrate(likelihood, prior) + julia> post isa MeasureBase.DensityMeasure + true + julia> logdensity_rel(post, Lebesgue(), 2) + -4.337877066409345 """ struct Likelihood{K,X} <: AbstractLikelihood k::K x::X - Likelihood(k::K, x::X) where {K<:AbstractTransitionKernel,X} = new{K,X}(k, x) - Likelihood(k::K, x::X) where {K<:Function,X} = new{K,X}(k, x) - Likelihood(μ, x) = Likelihood(kernel(μ), x) + Likelihood{K,X}(k, x) where {K,X} = new{K,X}(k, x) end +# For type stability, in case k is a type (resp. a constructor): +Likelihood(k, x::X) where {X} = Likelihood{Core.Typeof(k),X}(k, x) + (lik::AbstractLikelihood)(p) = exp(ULogarithmic, logdensityof(lik.k(p), lik.x)) DensityInterface.DensityKind(::AbstractLikelihood) = IsDensity() @@ -150,58 +127,86 @@ end export likelihoodof -""" - likelihoodof(k::AbstractTransitionKernel, x; constraints...) - likelihoodof(k::AbstractTransitionKernel, x, constraints::NamedTuple) +@doc raw""" + likelihoodof(k, x) -A likelihood is *not* a measure. Rather, a likelihood acts on a measure, through -the "pointwise product" `⊙`, yielding another measure. -""" -function likelihoodof end +Returns the likelihood of observing `x` under a family of probability +measures that is generated by a transition kernel `k(θ)`. + +`k(θ)` maps points in the parameter space to measures (resp. objects that can +be converted to measures) on a implicit set `Χ` that contains values like `x`. -likelihoodof(k, x, ::NamedTuple{()}) = Likelihood(k, x) +`likelihoodof(k, x)` returns a likelihood object. A likelihhood is **not** a +measure, it is a function from the parameter space to `ℝ₊`. Likelihood +objects can also be interpreted as "generic densities" (but **not** as +probability densities). -likelihoodof(k, x; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs)) +`likelihoodof(k, x)` implicitly chooses `ξ = rootmeasure(k(θ))` as the +reference measure on the observation set `Χ`. Note that this implicit +`ξ` **must** be independent of `θ`. -likelihoodof(k, x, pars::NamedTuple) = likelihoodof(kernel(k, pars), x) +`ℒₓ = likelihoodof(k, x)` has the mathematical interpretation -likelihoodof(k::AbstractTransitionKernel, x) = Likelihood(k, x) +```math +\mathcal{L}_x(\theta) = \frac{\rm{d}\, k(\theta)}{\rm{d}\, \chi}(x) +``` -export log_likelihood_ratio +`likelihoodof` must return an object that implements the +[`DensityInterface`](https://github.com/JuliaMath/DensityInterface.jl)` API +and `ℒₓ = likelihoodof(k, x)` must satisfy +```julia +log(ℒₓ(θ)) == logdensityof(ℒₓ, θ) ≈ logdensityof(k(θ), x) + +DensityKind(ℒₓ) isa IsDensity +``` + +By default, an instance of [`MeasureBase.Likelihood`](@ref) is returned. """ - log_likelihood_ratio(ℓ::Likelihood, p, q) +function likelihoodof end -Compute the log of the likelihood ratio, in order to compare two choices for -parameters. This is computed as +likelihoodof(k, x) = Likelihood(k, x) - logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x) +############################################################################### +# At the least, we need to think through in some more detail whether +# (log-)likelihood ratios expressed in this way are correct and useful. For now +# this code is commented out; we may remove it entirely in the future. -Since `logdensity_rel` can leave common base measure unevaluated, this can be -more efficient than +# export log_likelihood_ratio - logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x) -""" -log_likelihood_ratio(ℓ::Likelihood, p, q) = logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x) +# """ +# log_likelihood_ratio(ℓ::Likelihood, p, q) -# likelihoodof(k, x; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs)) +# Compute the log of the likelihood ratio, in order to compare two choices for +# parameters. This is computed as -export likelihood_ratio +# logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x) -""" - likelihood_ratio(ℓ::Likelihood, p, q) +# Since `logdensity_rel` can leave common base measure unevaluated, this can be +# more efficient than -Compute the log of the likelihood ratio, in order to compare two choices for -parameters. This is equal to +# logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x) +# """ +# log_likelihood_ratio(ℓ::Likelihood, p, q) = logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x) - density_rel(ℓ.k(p), ℓ.k(q), ℓ.x) +# # likelihoodof(k, x; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs)) -but is computed using LogarithmicNumbers.jl to avoid underflow and overflow. -Since `density_rel` can leave common base measure unevaluated, this can be -more efficient than +# export likelihood_ratio - logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x) -""" -function likelihood_ratio(ℓ::Likelihood, p, q) - exp(ULogarithmic, logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)) -end +# """ +# likelihood_ratio(ℓ::Likelihood, p, q) + +# Compute the log of the likelihood ratio, in order to compare two choices for +# parameters. This is equal to + +# density_rel(ℓ.k(p), ℓ.k(q), ℓ.x) + +# but is computed using LogarithmicNumbers.jl to avoid underflow and overflow. +# Since `density_rel` can leave common base measure unevaluated, this can be +# more efficient than + +# logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x) +# """ +# function likelihood_ratio(ℓ::Likelihood, p, q) +# exp(ULogarithmic, logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)) +# end diff --git a/src/combinators/pointwise.jl b/src/combinators/pointwise.jl deleted file mode 100644 index 778e7f4e..00000000 --- a/src/combinators/pointwise.jl +++ /dev/null @@ -1,30 +0,0 @@ -export ⊙ - -struct PointwiseProductMeasure{P,L} <: AbstractMeasure - prior::P - likelihood::L -end - -iterate(p::PointwiseProductMeasure, i = 1) = iterate((p.prior, p.likelihood), i) - -function Pretty.tile(d::PointwiseProductMeasure) - Pretty.pair_layout(Pretty.tile(d.prior), Pretty.tile(d.likelihood), sep = " ⊙ ") -end - -⊙(prior, ℓ) = pointwiseproduct(prior, ℓ) - -@inbounds function insupport(d::PointwiseProductMeasure, p) - prior, ℓ = d - istrue(insupport(prior, p)) && istrue(insupport(ℓ, p)) -end - -@inline function logdensity_def(d::PointwiseProductMeasure, p) - prior, ℓ = d - unsafe_logdensityof(ℓ, p) -end - -basemeasure(d::PointwiseProductMeasure) = d.prior - -function gentype(d::PointwiseProductMeasure) - gentype(d.prior) -end diff --git a/src/combinators/powerweighted.jl b/src/combinators/powerweighted.jl deleted file mode 100644 index 47f50da4..00000000 --- a/src/combinators/powerweighted.jl +++ /dev/null @@ -1,37 +0,0 @@ -export ↑ - -struct PowerWeightedMeasure{M,A} <: AbstractMeasure - parent::M - exponent::A -end - -logdensity_def(d::PowerWeightedMeasure, x) = d.exponent * logdensity_def(d.parent, x) - -basemeasure(d::PowerWeightedMeasure, x) = basemeasure(d.parent, x)↑d.exponent - -basemeasure(d::PowerWeightedMeasure) = basemeasure(d.parent)↑d.exponent - -function powerweightedmeasure(d, α) - isone(α) && return d - PowerWeightedMeasure(d, α) -end - -(d::AbstractMeasure)↑α = powerweightedmeasure(d, α) - -insupport(d::PowerWeightedMeasure, x) = insupport(d.parent, x) - -function Base.show(io::IO, d::PowerWeightedMeasure) - print(io, d.parent, " ↑ ", d.exponent) -end - -function powerweightedmeasure(d::PowerWeightedMeasure, α) - powerweightedmeasure(d.parent, α * d.exponent) -end - -function powerweightedmeasure(d::WeightedMeasure, α) - weightedmeasure(α * d.logweight, powerweightedmeasure(d.base, α)) -end - -function Pretty.tile(d::PowerWeightedMeasure) - Pretty.pair_layout(Pretty.tile(d.parent), Pretty.tile(d.exponent), sep = " ↑ ") -end diff --git a/src/combinators/product.jl b/src/combinators/product.jl index cb7a0aaf..516678f5 100644 --- a/src/combinators/product.jl +++ b/src/combinators/product.jl @@ -167,19 +167,6 @@ function testvalue(::Type{T}, d::AbstractProductMeasure) where {T} _map(m -> testvalue(T, m), marginals(d)) end -export ⊗ - -""" - ⊗(μs::AbstractMeasure...) - -`⊗` is a binary operator for building product measures. This satisfies the law - -``` - basemeasure(μ ⊗ ν) == basemeasure(μ) ⊗ basemeasure(ν) -``` -""" -⊗(μs::AbstractMeasure...) = productmeasure(μs) - ############################################################################### # I <: Base.Generator diff --git a/src/combinators/transformedmeasure.jl b/src/combinators/transformedmeasure.jl index 803b404b..dab76d5f 100644 --- a/src/combinators/transformedmeasure.jl +++ b/src/combinators/transformedmeasure.jl @@ -140,7 +140,7 @@ end # pullback """ - pullback(f, μ, volcorr = WithVolCorr()) + pullbck(f, μ, volcorr = WithVolCorr()) A _pullback_ is a dual concept to a _pushforward_. While a pushforward needs a map _from_ the support of a measure, a pullback requires a map _into_ the @@ -152,8 +152,11 @@ in terms of the inverse function; the "forward" function is not used at all. In some cases, we may be focusing on log-density (and not, for example, sampling). To manually specify an inverse, call -`pullback(InverseFunctions.setinverse(f, finv), μ, volcorr)`. +`pullbck(InverseFunctions.setinverse(f, finv), μ, volcorr)`. """ -function pullback(f, μ, volcorr::TransformVolCorr = WithVolCorr()) - pushfwd(setinverse(inverse(f), f), μ, volcorr) +function pullbck(f, μ, volcorr::TransformVolCorr = WithVolCorr()) + PushforwardMeasure(inverse(f), f, μ, volcorr) end +export pullbck + +@deprecate pullback(f, μ, volcorr::TransformVolCorr = WithVolCorr()) pullbck(f, μ, volcorr) diff --git a/src/combinators/weighted.jl b/src/combinators/weighted.jl index db239b50..124662b6 100644 --- a/src/combinators/weighted.jl +++ b/src/combinators/weighted.jl @@ -46,9 +46,6 @@ end Base.:*(m::AbstractMeasure, k::Real) = k * m -≪(::M, ::WeightedMeasure{R,M}) where {R,M} = true -≪(::WeightedMeasure{R,M}, ::M) where {R,M} = true - gentype(μ::WeightedMeasure) = gentype(μ.base) insupport(μ::WeightedMeasure, x) = insupport(μ.base, x) diff --git a/src/density.jl b/src/density.jl index 4862dcb1..0f4c8e03 100644 --- a/src/density.jl +++ b/src/density.jl @@ -20,8 +20,7 @@ For measures `μ` and `ν`, `Density(μ,ν)` represents the _density function_ `dμ/dν`, also called the _Radom-Nikodym derivative_: https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem#Radon%E2%80%93Nikodym_derivative -Instead of calling this directly, users should call `density_rel(μ, ν)` or -its abbreviated form, `𝒹(μ,ν)`. +Instead of calling this directly, users should call `density_rel(μ, ν)`. """ struct Density{M,B} <: AbstractDensity μ::M @@ -32,16 +31,6 @@ Base.:∘(::typeof(log), d::Density) = logdensity_rel(d.μ, d.base) Base.log(d::Density) = log ∘ d -export 𝒹 - -""" - 𝒹(μ, base) - -Compute the density (Radom-Nikodym derivative) of μ with respect to `base`. This -is a shorthand form for `density_rel(μ, base)`. -""" -𝒹(μ, base) = density_rel(μ, base) - density_rel(μ, base) = Density(μ, base) (f::Density)(x) = density_rel(f.μ, f.base, x) @@ -73,16 +62,6 @@ Base.:∘(::typeof(exp), d::LogDensity) = density_rel(d.μ, d.base) Base.exp(d::LogDensity) = exp ∘ d -export log𝒹 - -""" - log𝒹(μ, base) - -Compute the log-density (Radom-Nikodym derivative) of μ with respect to `base`. -This is a shorthand form for `logdensity_rel(μ, base)` -""" -log𝒹(μ, base) = logdensity_rel(μ, base) - logdensity_rel(μ, base) = LogDensity(μ, base) (f::LogDensity)(x) = logdensity_rel(f.μ, f.base, x) @@ -98,12 +77,13 @@ DensityInterface.funcdensity(d::LogDensity) = throw(MethodError(funcdensity, (d, base :: B end -A `DensityMeasure` is a measure defined by a density or log-density with respect -to some other "base" measure. +A `DensityMeasure` is a measure defined by a density or log-density with +respect to some other "base" measure. -Users should not call `DensityMeasure` directly, but should instead call `∫(f, -base)` (if `f` is a density function or `DensityInterface.IsDensity` object) or -`∫exp(f, base)` (if `f` is a log-density function). +Users should not instantiate `DensityMeasure` directly, but should instead +call `mintegral_exp(f, base)` (if `f` is a density function or +`DensityInterface.IsDensity` object) or `mintegral_exp(f, base)` (if `f` +is a log-density function). """ struct DensityMeasure{F,B} <: AbstractMeasure f::F @@ -120,56 +100,84 @@ end end function Pretty.tile(μ::DensityMeasure{F,B}) where {F,B} - result = Pretty.literal("DensityMeasure ∫(") + result = Pretty.literal("mintegrate(") result *= Pretty.pair_layout(Pretty.tile(μ.f), Pretty.tile(μ.base); sep = ", ") result *= Pretty.literal(")") end -export ∫ +basemeasure(μ::DensityMeasure) = μ.base -""" - ∫(f, base::AbstractMeasure) +logdensity_def(μ::DensityMeasure, x) = logdensityof(μ.f, x) -Define a new measure in terms of a density `f` over some measure `base`. -""" -∫(f, base) = _densitymeasure(f, base, DensityKind(f)) +density_def(μ::DensityMeasure, x) = densityof(μ.f, x) -_densitymeasure(f, base, ::IsDensity) = DensityMeasure(f, base) -function _densitymeasure(f, base, ::HasDensity) - @error "`∫(f, base)` requires `DensityKind(f)` to be `IsDensity()` or `NoDensity()`." -end -_densitymeasure(f, base, ::NoDensity) = DensityMeasure(funcdensity(f), base) +@doc raw""" + mintegrate(f, μ::AbstractMeasure)::AbstractMeasure -export ∫exp +Returns a new measure that represents the indefinite +[integral](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem) +of `f` with respect to `μ`. -""" - ∫exp(f, base::AbstractMeasure) +`ν = mintegrate(f, μ)` generates a measure `ν` that has the mathematical +interpretation -Define a new measure in terms of a log-density `f` over some measure `base`. +math``` +\nu(A) = \int_A f(a) \, \rm{d}\mu(a) +``` """ -∫exp(f, base) = _logdensitymeasure(f, base, DensityKind(f)) +function mintegrate end +export mintegrate -function _logdensitymeasure(f, base, ::IsDensity) - @error "`∫exp(f, base)` is not valid when `DensityKind(f) == IsDensity()`. Use `∫(f, base)` instead." -end -function _logdensitymeasure(f, base, ::HasDensity) - @error "`∫exp(f, base)` is not valid when `DensityKind(f) == HasDensity()`." +mintegrate(f, μ::AbstractMeasure) = _mintegrate_impl(f, μ, DensityKind(f)) + +_mintegrate_impl(f, μ, ::IsDensity) = DensityMeasure(f, μ) +function _mintegrate_impl(f, μ, ::HasDensity) + throw( + ArgumentError( + "`mintegrate(f, mu)` requires `DensityKind(f)` to be `IsDensity()` or `NoDensity()`.", + ), + ) end -_logdensitymeasure(f, base, ::NoDensity) = DensityMeasure(logfuncdensity(f), base) +_mintegrate_impl(f, μ, ::NoDensity) = DensityMeasure(funcdensity(f), μ) -basemeasure(μ::DensityMeasure) = μ.base +@doc raw""" + mintegrate_exp(log_f, μ::AbstractMeasure) -logdensity_def(μ::DensityMeasure, x) = logdensityof(μ.f, x) +Given a function `log_f` that semantically represents the log of a function +`f`, `mintegrate` returns a new measure that represents the indefinite +[integral](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem) +of `f` with respect to `μ`. -density_def(μ::DensityMeasure, x) = densityof(μ.f, x) +`ν = mintegrate_exp(log_f, μ)` generates a measure `ν` that has the +mathematical interpretation -""" - rebase(μ, ν) - -Express `μ` in terms of a density over `ν`. Satisfies +math``` +\nu(A) = \int_A e^{log(f(a))} \, \rm{d}\mu(a) = \int_A f(a) \, \rm{d}\mu(a) ``` -basemeasure(rebase(μ, ν)) == ν -density(rebase(μ, ν)) == 𝒹(μ,ν) -``` + +Note that `exp(log_f(...))` is usually not run explicitly, calculations that +involve the resulting measure are typically performed in log-space, +internally. """ -rebase(μ, ν) = ∫(𝒹(μ, ν), ν) +function mintegrate_exp end +export mintegrate_exp + +function mintegrate_exp(log_f, μ::AbstractMeasure) + _mintegrate_exp_impl(log_f, μ, DensityKind(log_f)) +end + +function _mintegrate_exp_impl(log_f, μ, ::IsDensity) + throw( + ArgumentError( + "`mintegrate_exp(log_f, μ)` is not valid when `DensityKind(log_f) == IsDensity()`. Use `mintegral(log_f, μ)` instead.", + ), + ) +end +function _mintegrate_exp_impl(log_f, μ, ::HasDensity) + throw( + ArgumentError( + "`mintegrate_exp(log_f, μ)` is not valid when `DensityKind(log_f) == HasDensity()`.", + ), + ) +end +_mintegrate_exp_impl(log_f, μ, ::NoDensity) = DensityMeasure(logfuncdensity(log_f), μ) diff --git a/src/measure_operators.jl b/src/measure_operators.jl new file mode 100644 index 00000000..5822d4de --- /dev/null +++ b/src/measure_operators.jl @@ -0,0 +1,131 @@ +""" + module MeasureOperators + +Defines the following operators for measures: + +* `f ⋄ μ == pushfwd(f, μ)` + +* `μ ⊙ f == inverse(f) ⋄ μ` +""" +module MeasureOperators + +using MeasureBase: AbstractMeasure +using MeasureBase: pushfwd, pullbck, mbind, productmeasure +using MeasureBase: mintegrate, mintegrate_exp, density_rel, logdensity_rel +using InverseFunctions: inverse +using Reexport: @reexport + +@doc raw""" + ⋄(f, μ::AbstractMeasure) = pushfwd(f, μ) + +The `\\diamond` operator denotes a pushforward operation: `ν = f ⋄ μ` +generates a +[pushforward measure](https://en.wikipedia.org/wiki/Pushforward_measure). + +A common mathematical notation for a pushforward is ``f_*μ``, but as +there is no "subscript-star" operator in Julia, we use `⋄`. + +See [`pushfwd(f, μ)`](@ref) for details. + +Also see [`ν ⊙ f`](@ref), the pullback operator. +""" +⋄(f, μ::AbstractMeasure) = pushfwd(f, μ) +export ⋄ + +@doc raw""" + ⊙(ν::AbstractMeasure, f) = pullbck(f, ν) + +The `\\odot` operator denotes a pullback operation. + +See also [`pullbck(ν, f)`](@ref) for details. Note that `pullbck` takes it's +arguments in different order, in keeping with the Julia convention of +passing functions as the first argument. A pullback is mathematically the +precomposition of a measure `μ`` with the function `f` applied to sets. so +`⊙` takes the measure as the first and the function as the second argument, +as common in mathematical notation for precomposition. + +A common mathematical notation for pullback in measure theory is +``f \circ μ``, but as `∘` is used for function composition in Julia and as +`f` semantically acts point-wise on sets, we use `⊙`. + +Also see [f ⋄ μ](@ref), the pushforward operator. +""" +⊙(ν::AbstractMeasure, f) = pullbck(f, ν) +export ⊙ + +""" + μ ▷ k = mbind(k, μ) + +The `\\triangleright` operator denotes a measure monadic bind operation. + +A common operator choice for a monadics bind operator is `>>=` (e.g. in +the Haskell programming language), but this has a different meaning in +Julia and there is no close equivalent, so we use `▷`. + +See [`mbind(k, μ)`](@ref) for details. Note that `mbind` takes its +arguments in different order, in keeping with the Julia convention of +passing functions as the first argument. `▷`, on the other hand, takes +its arguments in the order common for monadic binds in functional +programming (like the Haskell `>>=` operator) and mathematics. +""" +▷(μ::AbstractMeasure, k) = mbind(k, μ) +export ▷ + +# ToDo: Use `⨂` instead of `⊗` for better readability? +""" + ⊗(μs::AbstractMeasure...) = productmeasure(μs) + +`⊗` is an operator for building product measures. + +See [`productmeasure(μs)`](@ref) for details. +""" +⊗(μs::AbstractMeasure...) = productmeasure(μs) +export ⊗ + +""" + ∫(f, μ::AbstractMeasure) = mintegrate(f, μ) + +Denotes an indefinite integral of the function `f` with respect to the +measure `μ`. + +See [`mintegrate(f, μ)`](@ref) for details. +""" +∫(f, μ::AbstractMeasure) = mintegrate(f, μ) +export ∫ + +""" + ∫exp(f, μ::AbstractMeasure) = mintegrate_exp(f, μ) + +Generates a new measure that is the indefinite integral of `exp` of `f` +with respect to the measure `μ`. + +See [`mintegrate_exp(f, μ)`](@ref) for details. +""" +∫exp(f, μ::AbstractMeasure) = mintegrate_exp(f, μ) +export ∫exp + +""" + 𝒹(ν, μ) = density_rel(ν, μ) + +Compute the density, i.e. the +[Radom-Nikodym derivative](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem) +of `ν`` with respect to `μ`. + +For details, see [`density_rel(ν, μ)`}(@ref). +""" +𝒹(ν, μ::AbstractMeasure) = density_rel(ν, μ) +export 𝒹 + +""" + log𝒹(ν, μ) = logdensity_rel(ν, μ) + +Compute the log-density, i.e. the logarithm of the +[Radom-Nikodym derivative](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem) +of `ν`` with respect to `μ`. + +For details, see [`logdensity_rel(ν, μ)`}(@ref). +""" +log𝒹(ν, μ::AbstractMeasure) = logdensity_rel(ν, μ) +export log𝒹 + +end # module MeasureOperators diff --git a/src/parameterized.jl b/src/parameterized.jl index 78e43995..8b1c8c88 100644 --- a/src/parameterized.jl +++ b/src/parameterized.jl @@ -127,14 +127,3 @@ params(::Type{PM}) where {N,PM<:ParameterizedMeasure{N}} = N function paramnames(μ, constraints::NamedTuple{N}) where {N} tuple((k for k in paramnames(μ) if k ∉ N)...) end - -############################################################################### -# kernelfactor - -function kernelfactor(::Type{P}) where {N,P<:ParameterizedMeasure{N}} - (constructorof(P), N) -end - -function kernelfactor(::P) where {N,P<:ParameterizedMeasure{N}} - (constructorof(P), N) -end diff --git a/src/static.jl b/src/static.jl index 81c528ba..8ec313fb 100644 --- a/src/static.jl +++ b/src/static.jl @@ -49,7 +49,9 @@ Returns the length of `x` as a dynamic or static integer. """ maybestatic_length(x) = length(x) maybestatic_length(x::AbstractUnitRange) = length(x) -function maybestatic_length(::Static.OptionallyStaticUnitRange{<:StaticInteger{A},<:StaticInteger{B}}) where {A,B} +function maybestatic_length( + ::Static.OptionallyStaticUnitRange{<:StaticInteger{A},<:StaticInteger{B}}, +) where {A,B} StaticInt{B - A + 1}() end diff --git a/test/measure_operators.jl b/test/measure_operators.jl new file mode 100644 index 00000000..a3adaa8f --- /dev/null +++ b/test/measure_operators.jl @@ -0,0 +1,24 @@ +using Test + +using MeasureBase: AbstractMeasure +using MeasureBase: StdExponential, StdLogistic, StdUniform +using MeasureBase: pushfwd, pullbck, mbind, productmeasure +using MeasureBase: mintegrate, mintegrate_exp, density_rel, logdensity_rel +using MeasureBase.MeasureOperators: ⋄, ⊙, ▷, ⊗, ∫, ∫exp, 𝒹, log𝒹 + +@testset "MeasureOperators" begin + μ = StdExponential() + ν = StdUniform() + k(σ) = pushfwd(x -> σ * x, StdNormal()) + μs = (StdExponential(), StdLogistic(), StdUniform()) + f = sqrt + + @test @inferred(f ⋄ μ) == pushfwd(f, μ) + @test @inferred(ν ⊙ f) == pullbck(f, ν) + @test @inferred(μ ▷ k) == mbind(k, μ) + @test @inferred(⊗(μs...)) == productmeasure(μs) + @test @inferred(∫(f, μ)) == mintegrate(f, μ) + @test @inferred(∫exp(f, μ)) == mintegrate_exp(f, μ) + @test @inferred(𝒹(ν, μ)) == density_rel(ν, μ) + @test @inferred(log𝒹(ν, μ)) == logdensity_rel(ν, μ) +end diff --git a/test/runtests.jl b/test/runtests.jl index 8a2deaa4..6cd33a41 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -22,7 +22,7 @@ include("static.jl") # # detect_ambiguities_options..., # ) -d = ∫exp(x -> -x^2, Lebesgue(ℝ)) +d = mintegrate_exp(x -> -x^2, Lebesgue(ℝ)) # function draw2(μ) # x = rand(μ) @@ -148,7 +148,7 @@ end @test logdensityof(Lebesgue()^3, 2) == logdensityof(Lebesgue()^(3, 1), (2, 0)) end -Normal() = ∫exp(x -> -0.5x^2, Lebesgue(ℝ)) +Normal() = mintegrate_exp(x -> -0.5x^2, Lebesgue(ℝ)) @testset "Half" begin HalfNormal() = Half(Normal()) @@ -159,12 +159,10 @@ end @testset "Likelihood" begin ℓ = Likelihood(3) do (μ,) - ∫exp(Lebesgue(ℝ)) do x + mintegrate_exp(Lebesgue(ℝ)) do x -(x - μ)^2 end end - - @inferred logdensityof(Lebesgue() ⊙ ℓ, 2.0) end # @testset "Likelihood" begin @@ -236,13 +234,13 @@ end @testset "Density measures and Radon-Nikodym" begin x = randn() f(x) = x^2 - @test log(𝒹(∫exp(f, Lebesgue()), Lebesgue())(x)) ≈ f(x) + @test log(density_rel(mintegrate_exp(f, Lebesgue()), Lebesgue())(x)) ≈ f(x) - let f = 𝒹(∫exp(x -> x^2, Lebesgue()), Lebesgue()) + let f = density_rel(mintegrate_exp(x -> x^2, Lebesgue()), Lebesgue()) @test log(f(x)) ≈ x^2 end - let f = log𝒹(∫exp(x -> x^2, Normal()), Normal()) + let f = logdensity_rel(mintegrate_exp(x -> x^2, Normal()), Normal()) @test f(x) ≈ x^2 end end @@ -253,3 +251,5 @@ include("smf.jl") include("combinators/weighted.jl") include("combinators/transformedmeasure.jl") + +include("measure_operators.jl") diff --git a/test/static.jl b/test/static.jl index a6c50db2..f618124b 100644 --- a/test/static.jl +++ b/test/static.jl @@ -11,7 +11,7 @@ import FillArrays @test static(2) isa MeasureBase.IntegerLike @test true isa MeasureBase.IntegerLike @test static(true) isa MeasureBase.IntegerLike - + @test @inferred(MeasureBase.one_to(7)) isa Base.OneTo @test @inferred(MeasureBase.one_to(7)) == 1:7 @test @inferred(MeasureBase.one_to(static(7))) isa Static.SOneTo @@ -19,10 +19,13 @@ import FillArrays @test @inferred(MeasureBase.fill_with(4.2, (7,))) == FillArrays.Fill(4.2, 7) @test @inferred(MeasureBase.fill_with(4.2, (static(7),))) == FillArrays.Fill(4.2, 7) - @test @inferred(MeasureBase.fill_with(4.2, (3, static(7)))) == FillArrays.Fill(4.2, 3, 7) + @test @inferred(MeasureBase.fill_with(4.2, (3, static(7)))) == + FillArrays.Fill(4.2, 3, 7) @test @inferred(MeasureBase.fill_with(4.2, (3:7,))) == FillArrays.Fill(4.2, (3:7,)) - @test @inferred(MeasureBase.fill_with(4.2, (static(3):static(7),))) == FillArrays.Fill(4.2, (3:7,)) - @test @inferred(MeasureBase.fill_with(4.2, (3:7, static(2):static(5)))) == FillArrays.Fill(4.2, (3:7, 2:5)) + @test @inferred(MeasureBase.fill_with(4.2, (static(3):static(7),))) == + FillArrays.Fill(4.2, (3:7,)) + @test @inferred(MeasureBase.fill_with(4.2, (3:7, static(2):static(5)))) == + FillArrays.Fill(4.2, (3:7, 2:5)) @test MeasureBase.maybestatic_length(MeasureBase.one_to(7)) isa Int @test MeasureBase.maybestatic_length(MeasureBase.one_to(7)) == 7