Dev (#79)

* scale.jl => weighted.jl

* scale.jl => weighted.jl

* some refactoring

* Poisson

* exp

* ForGenerator

* basemeasure(::ForArray)

* I may have intended to import Base.rand

* logdensity(d::ForGenerator, x)

* binomial

* density (Funny we didn't already have this)

* io = IOContext(io, :compact => true)

* include("exp.jl")

* rand methods for Binomial

* CountingMeasure bugfixes

* tests

* tests

* half distributions

* bugfix

* Binomial parameterized by probit_p

* Update Project.toml

* compat

* export WeightedMeasure

* probit

* compat

* loosen bounds on InfiniteArrays version

* loosen LazyArrays version bounds

* HalfNormal(σ) = HalfNormal(σ = σ)

* changing some dependencies around

* `distproxy` for delegating to Distributions and MCM

* distproxy

* Tullio

* deps

* roll back Tullio stuff for now

* drop LazyArrays

* version bump

* tullio

* logpdf(::AbstractMeasure)

* dependencies

* Make `distproxy` a little more flexible

* distproxy for student T

* elementwise product

* starting on likelihood

* temporary un-tullio

* drop Exp for now

* pushforward

* dirac logdensity

* updating binomial

* sampletype(::CountingMeasure)= Int

* likelihood

* transforms

* transforms

* transforms et al

* kernelize

* testvalue

* (::Type{$μ{N}})(nt::NamedTuple{N,T}) where {N,T} = $μ{N,T}(nt)

* refactoring `@neasure`

* more refactoring

* deps

* first shot at order statistics

* Likelihoods aren't necessarily measures

* oops left this out

* basemeasure(::Normal)

* fixing basemeasure for Pushforward and Pullback

* export testvalue

* LKJL

* bugfix

* regrouping

* some testvalue bugs

* update deps

* revert base measure for now

* keysort parameters

* too many underscores

* asparams

* empty tuple method for Parameterized measures

* paramsort

* keysort => paramsort

* @measure Binomial

* more normal methods

* tests

* spacing

* bugfix

* zero(float(x)) => zero(x)

* Update Project.toml

* deps

* Update Project.toml

Author: cscherrer

Project.toml

@@ -4,13 +4,15 @@ authors = ["Chad Scherrer <[email protected]> and contributors"]
 version = "0.5.0"
 
 [deps]
+Accessors = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 DistributionsAD = "ced4e74d-a319-5a8a-b0ac-84af2272839c"
+DynamicIterators = "6c76993d-992e-5bf1-9e63-34920a5a5a38"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
 IterTools = "c8e1da08-722c-5040-9ed9-7db0dc04731e"
-KeywordDispatch = "5888135b-5456-5c80-a1b6-c91ef8180460"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
@@ -23,25 +25,31 @@ RandomNumbers = "e6cf234a-135c-5ec9-84dd-332b85af5143"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 StatsFuns = "4c63d2b9-4356-54db-8cca-17b64c39e42c"
+TransformVariables = "84d833dd-6860-57f9-a1a7-6da5db126cff"
+Tullio = "bc48ee85-29a4-5162-ae0b-a64e1601d4bc"
 
 [compat]
+Accessors = "0.1"
 ArrayInterface = "2,3"
+ConcreteStructs = "0.2"
 Distributions = "0.23, 0.24"
 DistributionsAD = "0.6"
+DynamicIterators = "0.4"
 InfiniteArrays = "0.7, 0.8, 0.9, 0.10"
 IntervalSets = "0.5"
 IterTools = "1"
-KeywordDispatch = "0.3"
 MLStyle = "0.4"
 MacroTools = "0.5"
 MappedArrays = "0.3, 0.4"
 MonteCarloMeasurements = "0.10"
 NamedTupleTools = "0.13"
-NestedTuples = "0.1, 0.2, 0.3"
+NestedTuples = "0.3"
 RandomNumbers = "1"
 SpecialFunctions = "0.10, 1"
 StaticArrays = "0.12, 1"
 StatsFuns = "0.9"
+TransformVariables = "0.3"
+Tullio = "0.2"
 julia = "1.5"
 
 [extras]

scratch/orderstats.jl

@@ -0,0 +1,13 @@
+struct OrderStatistic{n,r,B} <: AbstractMeasure
+    base::B
+end
+
+function basemeasure(d::OrderStatistic{n,r})
+    WeightedMeasure(logbeta(n+1, r+1), d.base)
+end
+
+function logdensity(d::OrderStatistic{n,r}, x)
+    logF(x) = logcdf(d.base, x)
+    logFᶜ(x) = logccdf(d.base, x)
+    return (r - 1) * logF(x) + (n - r) * logFᶜ(x)
+end

src/MeasureTheory.jl

@@ -4,7 +4,7 @@ using Random
 
 using MLStyle
 using NestedTuples
-
+using TransformVariables
 import Base
 import Distributions
 const Dists = Distributions
@@ -13,13 +13,17 @@ export ≪
 export sampletype
 
 export AbstractMeasure
-using InfiniteArrays: ∞
+using InfiniteArrays
+
+const ∞ = InfiniteArrays.∞
 
 export ∞
 
 abstract type AbstractMeasure end
 
-sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))
+sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))
+
+# sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))
 
 """
     logdensity(μ::AbstractMeasure{X}, x::X)
@@ -32,6 +36,7 @@ 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")
@@ -44,6 +49,9 @@ include("combinators/superpose.jl")
 include("combinators/product.jl")
 include("combinators/for.jl")
 include("combinators/power.jl")
+include("combinators/elementwise.jl")
+include("combinators/transforms.jl")
+include("combinators/spikemixture.jl")
 include("distributions.jl")
 include("rand.jl")
 include("probability/dirac.jl")
@@ -60,8 +68,9 @@ include("probability/inverse-gamma.jl")
 include("probability/bernoulli.jl")
 include("probability/poisson.jl")
 include("probability/binomial.jl")
-include("combinators/spikemixture.jl")
+include("probability/LKJL.jl")
 include("density.jl")
+include("likelihood.jl")
 # include("pushforward.jl")
 include("kernel.jl")
 include("distproxy.jl")

src/basemeasures/counting.jl

@@ -17,9 +17,9 @@ isprimitive(::CountingMeasure) = true
 # sampletype(::CountingMeasure{ℝ₊}) = Float64
 # sampletype(::CountingMeasure{𝕀}) = Float64
 
-sampletype(::CountingMeasure{IntegerRange{lo,hi}}) where {lo, hi} = Int
-
+sampletype(::CountingMeasure) = Int
 
+testvalue(μ::CountingMeasure{X}) where {X} = testvalue(X)
 
 logdensity(::CountingMeasure, x) = zero(float(x))
 

src/basemeasures/lebesgue.jl

@@ -19,5 +19,8 @@ sampletype(::Lebesgue{ℝ}) = Float64
 sampletype(::Lebesgue{ℝ₊}) = Float64
 sampletype(::Lebesgue{𝕀}) = Float64
 
+testvalue(::Lebesgue{ℝ}) = 0.0
+testvalue(::Lebesgue{𝕀}) = 0.5
+testvalue(::Lebesgue{ℝ₊}) = 1.0
 
-logdensity(::Lebesgue, x) = zero(float(x))
+logdensity(::Lebesgue, x) = zero(x)

src/combinators/elementwise.jl

@@ -0,0 +1,58 @@
+export ⊙
+
+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 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))

src/combinators/product.jl

@@ -56,19 +56,54 @@ function Base.:*(μ::M, ν::N) where {M <: AbstractMeasure, N <: AbstractMeasure
     ProductMeasure(data...)
 end
 
-@propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure, x)
+using Tullio
+
+function logdensity(d::ProductMeasure{A}, x) where {T, A<:AbstractArray}
     data = d.data
     @boundscheck size(data) == size(x) || throw(BoundsError)
-
     s = 0.0
-    Δs(j) = @inbounds logdensity(data[j], x[j])
-
-    @inbounds @simd for j in eachindex(x)
-        s += Δs(j)
+    Δs(i) = logdensity(data[i], x[i])
+    for i in eachindex(x)
+        s += Δs(i)
     end
-    s
+    return s
 end
 
+# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,1}}
+#     data = d.data
+#     @boundscheck size(data) == size(x) || throw(BoundsError)
+#     @tullio s = logdensity(data[i], x[i])
+#     s
+# end
+
+# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,2}}
+#     data = d.data
+#     @boundscheck size(data) == size(x) || throw(BoundsError)
+#     @tullio s = @inbounds logdensity(data[i,j], x[i,j])
+#     s
+# end
+
+# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,3}}
+#     data = d.data
+#     @boundscheck size(data) == size(x) || throw(BoundsError)
+#     @tullio s = @inbounds logdensity(data[i,j,k], x[i,j,k])
+#     s
+# end
+
+# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,4}}
+#     data = d.data
+#     @boundscheck size(data) == size(x) || throw(BoundsError)
+#     @tullio s = @inbounds logdensity(data[i,j,k,l], x[i,j,k,l])
+#     s
+# end
+
+# @propagate_inbounds function MeasureTheory.logdensity(d::ProductMeasure{A}, x) where{T, A<:AbstractArray{T,5}}
+#     data = d.data
+#     @boundscheck size(data) == size(x) || throw(BoundsError)
+#     @tullio s = @inbounds logdensity(data[i,j,k,l,m], x[i,j,k,l,m])
+#     s
+# end
+
 export rand!
 using Random: rand!, GLOBAL_RNG, AbstractRNG
 

src/combinators/transforms.jl

@@ -64,7 +64,7 @@ TransformVariables.as(ν::Pushforward) = ν.f ∘ as(ν.μ)
 
 TransformVariables.as(μ::Pullback) = inverse(μ.f) ∘ μ.ν
 
-as(::Lebesgue) = asℝ
+TransformVariables.as(::Lebesgue) = asℝ
 
 # t = as𝕀
 # μ = Normal()

src/density.jl

@@ -86,10 +86,12 @@ function logdensity(μ::AbstractMeasure, ν::AbstractMeasure, x)
     result += _logdensity(basemeasure(bμ), basemeasure(bν), x)
 end
 
-_logdensity(::Lebesgue{ℝ}, ::Lebesgue{ℝ}, x) = zero(float(x))
+_logdensity(::Lebesgue{ℝ}, ::Lebesgue{ℝ}, x) = zero(x)
 
 export density
 
-density(μ::AbstractMeasure, ν::AbstractMeasure, x) = Exp(logdensity(μ, ν, x))
+density(μ::AbstractMeasure, ν::AbstractMeasure, x) = exp(logdensity(μ, ν, x))
 
-density(μ::AbstractMeasure, x) = Exp(logdensity(μ, x))
+density(μ::AbstractMeasure, x) = exp(logdensity(μ, x))
+
+Δlogdensity(μ::AbstractMeasure, x) = logdensity(μ,x)

src/distproxy.jl

@@ -8,6 +8,7 @@ PROXIES = Dict(
         :mean
         :std
         :entropy
+        :cdf
         ],
     :MonteCarloMeasurements => [
         :Particles
@@ -19,7 +20,7 @@ for m in keys(PROXIES)
         @eval begin
             import $m: $f
             export $f
-            $m.$f(d::AbstractMeasure) = $m.$f(MeasureTheory.distproxy(d))
+            $m.$f(d::AbstractMeasure, args...) = $m.$f(MeasureTheory.distproxy(d), args...)
         end
     end
 end

src/domains.jl

@@ -55,6 +55,8 @@ function Base.show(io::IO, r::IntegerRange{lo, hi}) where {lo, hi}
     print(io, "ℤ[", lo, ":", hi, "]")
 end
 
+testvalue(::IntegerRange{lo, hi}) where {lo, hi} = lo
+
 ###########################################################
 # Real intervals
 

src/kernel.jl

@@ -17,9 +17,18 @@ If the argument is a named tuple `(;a=f1, b=f1)`, `κ(x)` is defined as
 struct Kernel{T,S}
     ops::S
 end
+
+export kernel
 kernel(op, ::Type{M}) where {M} = Kernel{M,typeof(op)}(op)
 kernel(::Type{M}; ops...) where {M} = Kernel{M,typeof(ops.data)}(ops.data)
 mapcall(t, x) = map(func -> func(x), t)
 (k::Kernel{M,<:Tuple})(x) where {M} = M(mapcall(k.ops, x)...)
 (k::Kernel{M,<:NamedTuple})(x) where {M} = M(;mapcall(k.ops, x)...)
 (k::Kernel{M})(x) where {M} = M(k.ops(x)...)
+
+export kernelize
+
+function kernelize(μ::M) where {N, T, M <: ParameterizedMeasure{N,T}}
+    C = M.name.wrapper
+    (kernel(NamedTuple{N}, C), values(getfield(μ, :par)))
+end

src/likelihood.jl

@@ -0,0 +1,13 @@
+export Likelihood
+
+struct Likelihood{T,X}
+    x::X
+end
+
+Likelihood(T::Type, x::X) where {X} = Likelihood{T,X}(x)
+
+Likelihood(μ::T, x::X) where {X, T<:AbstractMeasure} = Likelihood{T,X}(x)
+
+function logdensity(μ::Likelihood{T,X}, p) where {T,X}
+    logdensity(T(p), μ.x)
+end

src/macros.jl

@@ -41,32 +41,21 @@ macro capture(template, ex, action)
 end
 
 function _measure(expr)
-    @capture $rel($μ($(p...)), $base) expr begin
-        base = replace(s -> s∈p, s -> :(μ.$s), base)
-
+    @capture $μ($(p...)) expr begin
         q = quote
             struct $μ{N,T} <: ParameterizedMeasure{N,T}
                 par :: NamedTuple{N,T}
             end
 
-            function MeasureTheory.basemeasure(μ::$μ{P}) where {P}
-                return $base
-            end
- 
-            (::$μ{P} ≪ ::typeof(representative($base))) where {P}  = true
-        end    
-
-
+            (::Type{$μ{N}})(nt::NamedTuple{N,T}) where {N,T} = $μ{N,T}(nt) 
+        end   
+        
         if !isempty(p)
             # e.g. Normal(μ,σ) = Normal(;μ=μ, σ=σ)
             # Requires Julia 1.5
             push!(q.args, :($μ($(p...)) = $μ(;$(p...))))
         end
-    
-        if rel == (:≃)
-            push!(q.args, :((::typeof(representative($base)) ≪ ::$μ{P}) where {P} = true))
-        end
-    
+        
         return q
     end
 end