Introduce mintegrate and mintegrate_exp

Removes the integral operators from MeasureBase,
to be re-introduced in the submodule MeasureOperators.

Also improves the likelihood documentation.

Author: oschulz

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
@@ -91,34 +91,35 @@ Similarly to the above, we have
 
 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 actually the Lebesgue integral,
+of the likelihood with respect to the prior.
 
 For example, say we have
 
     μ ~ Normal()
     x ~ Normal(μ,σ)
     σ = 1
 
-and we observe `x=3`. We can compute the posterior measure on `μ` as
-
-    julia> post = Normal() ⊙ Likelihood(Normal{(:μ, :σ)}, (σ=1,), 3)
-    Normal() ⊙ Likelihood(Normal{(:μ, :σ), T} where T, (σ = 1,), 3)
+and we observe `x=3`. We can compute the (non-normalized) posterior measure on
+`μ` as
 
-    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::K, x::X) where {K,X} = new{K,X}(k, x)
+#!!!!!!!!!!!    # For type stability if `K isa UnionAll (e.g. a parameterized MeasureType)`
+    Likelihood(::Type{K}, x::X) where {K<:AbstractMeasure,X} = new{K,X}(K, x)
 end
 
 (lik::AbstractLikelihood)(p) = exp(ULogarithmic, logdensityof(lik.k(p), lik.x))
@@ -150,58 +151,87 @@ 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)` 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, ::NamedTuple{()}) = Likelihood(k, x)
+`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; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs))
+`ℒₓ = likelihoodof(k, x)` has the mathematical interpretation
 
-likelihoodof(k, x, pars::NamedTuple) = likelihoodof(kernel(k, pars), x)
+```math
+\mathcal{L}_x(\theta) = \frac{\rm{d}\, k(\theta)}{\rm{d}\, \chi}(x)
+```
 
-likelihoodof(k::AbstractTransitionKernel, x) = Likelihood(k, x)
+`likelihoodof` must return an object that implements the
+[`DensityInterface`](https://github.com/JuliaMath/DensityInterface.jl)` API
+and `ℒₓ = likelihoodof(k, x)` must satisfy
 
-export log_likelihood_ratio
+```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)
 
-Since `logdensity_rel` can leave common base measure unevaluated, this can be
-more efficient than
+###############################################################################
+# 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.
 
-    logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x)
-"""
-log_likelihood_ratio(ℓ::Likelihood, p, q) = logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)
+# export log_likelihood_ratio
 
-# likelihoodof(k, x; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs))
+# """
+#     log_likelihood_ratio(ℓ::Likelihood, p, q)
 
-export likelihood_ratio
+# Compute the log of the likelihood ratio, in order to compare two choices for
+# parameters. This is computed as
 
-"""
-    likelihood_ratio(ℓ::Likelihood, p, q)
+#     logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x)
 
-Compute the log of the likelihood ratio, in order to compare two choices for
-parameters. This is equal to
+# Since `logdensity_rel` can leave common base measure unevaluated, this can be
+# more efficient than
 
-    density_rel(ℓ.k(p), ℓ.k(q), ℓ.x)
+#     logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x)
+# """
+# log_likelihood_ratio(ℓ::Likelihood, p, q) = logdensity_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
+# # likelihoodof(k, x; kwargs...) = likelihoodof(k, x, NamedTuple(kwargs))
 
-    logdensityof(ℓ.k(p), ℓ.x) - logdensityof(ℓ.k(q), ℓ.x)
-"""
-function likelihood_ratio(ℓ::Likelihood, p, q)
-    exp(ULogarithmic, logdensity_rel(ℓ.k(p), ℓ.k(q), ℓ.x))
-end
+# export likelihood_ratio
+
+# """
+#     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

src/density.jl

@@ -98,12 +98,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,48 +121,77 @@ 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
+
+logdensity_def(μ::DensityMeasure, x) = logdensityof(μ.f, x)
+
+density_def(μ::DensityMeasure, x) = densityof(μ.f, x)
 
-"""
-    ∫(f, base::AbstractMeasure)
 
-Define a new measure in terms of a density `f` over some measure `base`.
+
+@doc raw"""
+    mintegrate(f, μ::AbstractMeasure)::AbstractMeasure
+
+Returns a new measure that represents the indefinite
+[integral](https://en.wikipedia.org/wiki/Radon%E2%80%93Nikodym_theorem)
+of `f` with respect to `μ`.
+
+`ν = mintegrate(f, μ)` generates a measure `ν` that has the mathematical
+interpretation
+
+math```
+\nu(A) = \int_A f(a) \, \rm{d}\mu(a)
+```
 """
-∫(f, base) = _densitymeasure(f, base, DensityKind(f))
+function mintegrate end
+export mintegrate
+
+mintegrate(f, μ::AbstractMeasure) = _mintegrate_impl(f, μ, DensityKind(f))
 
-_densitymeasure(f, base, ::IsDensity) = DensityMeasure(f, base)
-function _densitymeasure(f, base, ::HasDensity)
-    @error "`∫(f, base)` requires `DensityKind(f)` to be `IsDensity()` or `NoDensity()`."
+_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
-_densitymeasure(f, base, ::NoDensity) = DensityMeasure(funcdensity(f), base)
+_mintegrate_impl(f, μ, ::NoDensity) = DensityMeasure(funcdensity(f), μ)
 
-export ∫exp
 
-"""
-    ∫exp(f, base::AbstractMeasure)
+@doc raw"""
+    mintegrate_exp(log_f, μ::AbstractMeasure)
+
+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 `μ`.
+
+`ν = mintegrate_exp(log_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 e^{log(f(a))} \, \rm{d}\mu(a) = \int_A f(a) \, \rm{d}\mu(a)
+```
+
+Note that `exp(log_f(...))` is usually not run explicitly, calculations that
+involve the resulting measure are typically performed in log-space,
+internally.
 """
-∫exp(f, base) = _logdensitymeasure(f, base, DensityKind(f))
+function mintegrate_exp end
+export mintegrate_exp
+
+mintegrate_exp(log_f, μ::AbstractMeasure) = _mintegrate_exp_impl(log_f, μ, DensityKind(log_f))
 
-function _logdensitymeasure(f, base, ::IsDensity)
-    @error "`∫exp(f, base)` is not valid when `DensityKind(f) == IsDensity()`. Use `∫(f, base)` instead."
+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 _logdensitymeasure(f, base, ::HasDensity)
-    @error "`∫exp(f, base)` is not valid when `DensityKind(f) == HasDensity()`."
+function _mintegrate_exp_impl(log_f, μ, ::HasDensity)
+    throw(ArgumentError("`mintegrate_exp(log_f, μ)` is not valid when `DensityKind(log_f) == HasDensity()`."))
 end
-_logdensitymeasure(f, base, ::NoDensity) = DensityMeasure(logfuncdensity(f), base)
+_mintegrate_exp_impl(log_f, μ, ::NoDensity) = DensityMeasure(logfuncdensity(log_f), μ)
 
-basemeasure(μ::DensityMeasure) = μ.base
-
-logdensity_def(μ::DensityMeasure, x) = logdensityof(μ.f, x)
-
-density_def(μ::DensityMeasure, x) = densityof(μ.f, x)
 
 """
     rebase(μ, ν)
@@ -172,4 +202,4 @@ basemeasure(rebase(μ, ν)) == ν
 density(rebase(μ, ν)) == 𝒹(μ,ν)
 ``` 
 """
-rebase(μ, ν) = ∫(𝒹(μ, ν), ν)
+rebase(μ, ν) = mintegrate(density_rel(μ, ν), ν)

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,7 +159,7 @@ end
 
 @testset "Likelihood" begin
     ℓ = Likelihood(3) do (μ,)
-        ∫exp(Lebesgue(ℝ)) do x
+        mintegrate_exp(Lebesgue(ℝ)) do x
             -(x - μ)^2
         end
     end
@@ -236,13 +236,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(𝒹(mintegrate_exp(f, Lebesgue()), Lebesgue())(x)) ≈ f(x)
 
-    let f = 𝒹(∫exp(x -> x^2, Lebesgue()), Lebesgue())
+    let f = 𝒹(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 = log𝒹(mintegrate_exp(x -> x^2, Normal()), Normal())
         @test f(x) ≈ x^2
     end
 end