populate from MeasureTheory

Author: cscherrer

Project.toml

@@ -3,6 +3,15 @@ uuid = "fa1605e6-acd5-459c-a1e6-7e635759db14"
 authors = ["Chad Scherrer <[email protected]> and contributors"]
 version = "0.1.0"
 
+[deps]
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
+MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
+MappedArrays = "dbb5928d-eab1-5f90-85c2-b9b0edb7c900"
+NamedTupleTools = "d9ec5142-1e00-5aa0-9d6a-321866360f50"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+
 [compat]
 julia = "1"
 

src/MeasureBase.jl

@@ -1,5 +1,68 @@
 module MeasureBase
 
-# Write your package code here.
+
+using Random
+
+using ConcreteStructs
+using MLStyle
+
+export ≪
+export sampletype
+
+export AbstractMeasure
+
+abstract type AbstractMeasure end
+
+sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))
+
+# sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))
+
+"""
+    logdensity(μ::AbstractMeasure{X}, x::X)
+
+Compute the logdensity of the measure μ at the point x. This is the standard way
+to define `logdensity` for a new measure. the base measure is implicit here, and
+is understood to be `basemeasure(μ)`.
+
+Methods for computing density relative to other measures will be
+"""
+function logdensity end
+
+include("paramorder.jl")
+include("exp.jl")
+include("domains.jl")
+include("utils.jl")
+include("absolutecontinuity.jl")
+include("basemeasures.jl")
+include("parameterized.jl")
+include("macros.jl")
+include("combinators/weighted.jl")
+include("combinators/superpose.jl")
+include("combinators/product.jl")
+include("combinators/for.jl")
+include("combinators/power.jl")
+include("combinators/likelihood.jl")
+include("combinators/elementwise.jl")
+include("combinators/spikemixture.jl")
+include("rand.jl")
+include("probability/dirac.jl")
+include("probability/normal.jl")
+include("probability/studentt.jl")
+include("probability/cauchy.jl")
+include("probability/laplace.jl")
+include("probability/uniform.jl")
+include("probability/beta.jl")
+include("probability/gumbel.jl")
+include("probability/exponential.jl")
+# include("probability/mvnormal.jl")
+include("probability/inverse-gamma.jl")
+include("probability/bernoulli.jl")
+include("probability/poisson.jl")
+include("probability/binomial.jl")
+# include("probability/LKJL.jl")
+include("probability/negativebinomial.jl")
+include("density.jl")
+# include("pushforward.jl")
+include("kernel.jl")
 
 end

src/MeasureTheory.jl

@@ -0,0 +1,3 @@
+module MeasureTheory
+
+end # module

src/absolutecontinuity.jl

@@ -0,0 +1,78 @@
+
+"""
+    ≪(μ,ν)
+
+# Absolute continuity
+
+A measure μ is _absolutely continuous_ with respect to ν, written μ ≪ ν, if
+ν(A)==0 implies μ(A)==0 for every ν-measurable set A.
+
+Less formally, suppose we have a set A with ν(A)==0. If μ(A)≠0, then there can
+be no way to "reweight" ν to get to μ. We can't make something from nothing.
+
+This "reweighting" is really a density function. If μ≪ν, then there is some
+function f that makes `μ == ∫(f,ν)` (see the help section for `∫`).
+
+We can get this f directly via the Radon-Nikodym derivative, `f == 𝒹(μ,ν)` (see
+the help section for `𝒹`).
+
+Note that `≪` is not a partial order, because it is not antisymmetric. That is
+to say, it's possible (in fact, common) to have two different measures `μ` and
+`ν` with `μ ≪ ν` and `ν ≪ μ`. A simple example of this is 
+```
+μ = Normal()
+ν = Lebesgue(ℝ)
+```
+
+When `≪` holds in both directions, the measures μ and ν are _equivalent_,
+written `μ ≃ ν`. See the help section for `≃` for more information.
+"""
+function ≪ end
+
+
+export ≃
+
+"""
+    ≃(μ,ν)
+
+# Equivalence of Measure
+
+Measures μ and ν on the same space X are equivalent, written `μ ≃ ν`, if `μ ≪ ν`
+and `ν ≪ μ`. Note that this is often written `~` in the literature, but this is
+overloaded in probabilistic programming, so we use this alternate notation. 
+
+Also note that equivalence is very different from equality. For two equivalent
+measures, the sets of non-zero measure will be identical, but what that measure
+is in each case can be very different. 
+"""
+function ≃(μ,ν)
+    return (μ≪ν && ν≪μ)
+end
+
+export representative
+
+"""
+    representative(μ::AbstractMeasure) -> AbstractMeasure
+"""
+function representative(μ)
+    function f(μ)
+        # Check if we're done
+        isprimitive(μ) && return μ
+        
+        ν = basemeasure(μ)
+
+        # TODO: Make sure we don't leave the equivalence class
+        # Make sure not to leave the equivalence class
+        # (ν ≪ μ) || return μ
+
+        return ν
+    end
+
+    fix(f, μ)
+end
+
+# TODO: ≪ needs more work
+function ≪(μ, ν)
+    μ == ν && return true
+    representative(μ) ≪ representative(ν) && return true
+end

src/basemeasures.jl

@@ -0,0 +1,35 @@
+export basemeasure
+
+"""
+    basemeasure(μ)
+
+Many measures are defined in terms of a logdensity relative to some base
+measure. This makes it important to be able to find that base measure.
+
+For measures not defined in this way, we'll typically have `basemeasure(μ) == μ`.
+"""
+function basemeasure end
+
+include("basemeasures/trivial.jl")
+include("basemeasures/lebesgue.jl")
+include("basemeasures/counting.jl")
+
+
+
+export isprimitive
+
+"""
+    isprimitive(μ)
+
+Most measures are defined in terms of other measures, for example using a
+density or a pushforward. Those that are not are considered (in this library,
+it's not a general measure theory thing) to be _primitive_. The canonical
+example of a primitive measure is `Lebesgue(X)` for some `X`.
+
+The default method is
+    isprimitive(μ) = false
+
+So when adding a new primitive measure, it's necessary to add a method for its type
+that returns `true`.
+"""
+isprimitive(μ) = false

src/basemeasures/counting.jl

@@ -0,0 +1,26 @@
+export CountingMeasure
+
+struct CountingMeasure{X} <: AbstractMeasure end
+
+function Base.show(io::IO, μ::CountingMeasure{X}) where {X}
+    io = IOContext(io, :compact => true)
+    print(io, "CountingMeasure(", X, ")")
+end
+
+CountingMeasure(X) = CountingMeasure{X}()
+
+basemeasure(μ::CountingMeasure) = μ
+
+isprimitive(::CountingMeasure) = true
+
+# sampletype(::CountingMeasure{ℝ}) = Float64
+# sampletype(::CountingMeasure{ℝ₊}) = Float64
+# sampletype(::CountingMeasure{𝕀}) = Float64
+
+sampletype(::CountingMeasure) = Int
+
+testvalue(μ::CountingMeasure{X}) where {X} = testvalue(X)
+
+logdensity(::CountingMeasure, x) = zero(float(x))
+
+# (::CountingMeaure)(s) = length(Set(s))

src/basemeasures/lebesgue.jl

@@ -0,0 +1,27 @@
+# Lebesgue measure
+
+export Lebesgue
+
+struct Lebesgue{X} <: AbstractMeasure end
+
+function Base.show(io::IO, μ::Lebesgue{X}) where X
+    io = IOContext(io, :compact => true)
+    print(io, "Lebesgue(", X, ")")
+end
+
+Lebesgue(X) = Lebesgue{X}()
+
+basemeasure(μ::Lebesgue) = μ
+
+isprimitive(::Lebesgue) = true
+
+sampletype(::Lebesgue{ℝ}) = Float64
+sampletype(::Lebesgue{ℝ₊}) = Float64
+sampletype(::Lebesgue{𝕀}) = Float64
+
+testvalue(::Lebesgue{ℝ}) = 0.0
+testvalue(::Lebesgue{𝕀}) = 0.5
+testvalue(::Lebesgue{ℝ₊}) = 1.0
+testvalue(::Lebesgue{<:Real}) = 0.0
+
+logdensity(::Lebesgue, x) = zero(x)

src/basemeasures/trivial.jl

@@ -0,0 +1,5 @@
+export TrivialMeasure
+
+struct TrivialMeasure <: AbstractMeasure end
+
+sampletype(::TrivialMeasure) = Nothing

src/combinators/elementwise.jl

@@ -0,0 +1,63 @@
+export ⊙
+
+@concrete terse struct ElementwiseProductMeasure{T} <: AbstractMeasure
+    data :: T
+
+    ElementwiseProductMeasure(μs...) = new{typeof(μs)}(μs)
+    ElementwiseProductMeasure(μs) = new{typeof(μs)}(μs)
+end
+
+Base.size(μ::ElementwiseProductMeasure) = size(μ.data)
+
+function Base.show(io::IO, μ::ElementwiseProductMeasure)
+    io = IOContext(io, :compact => true)
+    print(io, join(string.(μ.data), " ⊙ "))
+end
+
+function Base.show_unquoted(io::IO, μ::ElementwiseProductMeasure, indent::Int, prec::Int)
+    if Base.operator_precedence(:*) ≤ prec
+        print(io, "(")
+        show(io, μ)
+        print(io, ")")
+    else
+        show(io, μ)
+    end
+    return nothing
+end
+
+Base.length(m::ElementwiseProductMeasure{T}) where {T} = length(m.data)
+
+function ⊙(μ::ElementwiseProductMeasure{X}, ν::ElementwiseProductMeasure{Y}) where {X,Y}
+    data = (μ.data..., ν.data...)
+    ElementwiseProductMeasure(data...)
+end
+
+function ⊙(μ, ν::ElementwiseProductMeasure{Y}) where {Y}
+    data = (μ, ν.data...)
+    ElementwiseProductMeasure(data...)
+end
+
+function ⊙(μ::ElementwiseProductMeasure{X}, ν::N) where {X, N <: AbstractMeasure}
+    data = (μ.data..., ν)
+    ElementwiseProductMeasure(data...)
+end
+
+function ⊙(μ::M, ν::N) where {M <: AbstractMeasure, N <: AbstractMeasure}
+    data = (μ, ν)
+    ElementwiseProductMeasure(data...)
+end
+
+function ⊙(μ::AbstractMeasure, ℓ::LogLikelihood)
+    data = (μ, ℓ)
+    ElementwiseProductMeasure(data...)
+end
+
+function logdensity(d::ElementwiseProductMeasure, x)
+    sum((logdensity(dⱼ, x) for dⱼ in d.data))
+end
+
+function sampletype(d::ElementwiseProductMeasure) 
+    @inbounds sampletype(first(d.data))
+end
+
+basemeasure(μ::ElementwiseProductMeasure) =  @inbounds basemeasure(first(d.data))