Changes for MeasureBase v0.15

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

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

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

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,63 +64,40 @@ 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
 
     μ ~ Normal()
     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