analytical form for convolution and bind of two Normal distributions

src/MeasureTheory.jl

@@ -58,7 +58,8 @@ import MeasureBase:
     paramnames,
     ∫,
     𝒹,
-    ∫exp
+    ∫exp,
+    bind
 import MeasureBase: ≪
 using MeasureBase: BoundedInts, BoundedReals, CountingMeasure, IntegerDomain, IntegerNumbers
 using MeasureBase: weightedmeasure, restrict
@@ -127,6 +128,7 @@ include("parameterized.jl")
 
 include("macros.jl")
 include("combinators/affine.jl")
+include("combinators/convolve.jl")
 include("combinators/weighted.jl")
 include("combinators/product.jl")
 include("combinators/transforms.jl")
@@ -165,4 +167,7 @@ include("transforms/corrcholesky.jl")
 include("transforms/ordered.jl")
 
 include("distproxy.jl")
+
+include("parameterized/pairwise/normal_normal.jl")
+
 end # module

src/combinators/affine.jl

@@ -36,7 +36,7 @@ import InverseFunctions: inverse
 @inline inverse(f::AffineTransform{(:λ,)}) = AffineTransform((σ = f.λ,))
 @inline inverse(f::AffineTransform{(:μ,)}) = AffineTransform((μ = -f.μ,))
 
-# `size(f) == (m,n)` means `f : ℝⁿ → ℝᵐ`  
+# `size(f) == (m,n)` means `f : ℝⁿ → ℝᵐ`
 Base.size(f::AffineTransform{(:μ, :σ)}) = size(f.σ)
 Base.size(f::AffineTransform{(:μ, :λ)}) = size(f.λ)
 Base.size(f::AffineTransform{(:σ,)}) = size(f.σ)
@@ -339,3 +339,16 @@ end
 @inline function Distributions.cdf(d::Affine, x)
     cdf(parent(d), inverse(d.f)(x))
 end
+
+function mean(d::Affine)
+    m = mean(parent(d))
+    f = getfield(d, :f)
+    return f(m)
+end
+
+# std only for univariate distributions
+std(d::Affine{(:μ,)}) = std(parent(d))
+std(d::Affine{(:σ,)}) = d.σ * std(parent(d))
+std(d::Affine{(:λ,)}) = d.λ \ std(parent(d))
+std(d::Affine{(:μ, :σ)}) = d.σ * std(parent(d))
+std(d::Affine{(:μ, :λ)}) = d.λ \ std(parent(d))

src/combinators/convolve.jl

@@ -0,0 +1,16 @@
+struct Convolution{M,N} <: AbstractMeasure
+    μ::M
+    ν::N
+end
+
+"""
+If μ, ν are subtypes of `AbstractMeasure` or satisfy the Measure interface,
+then `convolve(μ, ν)` is a measure, called the convolution of μ and ν.
+"""
+convolve(μ, ν) = Convolution(μ, ν)
+
+function Base.rand(rng::AbstractRNG, ::Type{T}, d::Convolution) where {T}
+    x = rand(rng, T, d.μ)
+    y = rand(rng, T, d.ν)
+    return x+y
+end

src/parameterized/pairwise/normal_normal.jl

@@ -0,0 +1,13 @@
+function convolve(μ::Normal, ν::Normal)
+    Normal(mean(μ) + mean(ν), hypot(std(μ), std(ν)))
+end
+
+function bind(μ::Normal,
+    ::ParameterizedTransitionKernel{Type{Normal{(:μ,)}}, typeof(identity), (:μ,), Tuple{typeof(identity)}})
+   convolve(μ, Normal())
+end
+
+function bind(μ::Normal,
+    k::ParameterizedTransitionKernel{Type{Normal{(:μ, :σ)}}, typeof(identity), (:μ, :σ), Tuple{typeof(identity), T}} where T<:Number)
+   convolve(μ, Normal(σ=k.param_maps.σ))
+end

test/runtests.jl

@@ -277,7 +277,7 @@ end
 
 @testset "Product of Diracs" begin
     x = randn(3)
-    t = as(productmeasure(Dirac.(x))) 
+    t = as(productmeasure(Dirac.(x)))
     @test transform(t, []) == x
 end
 
@@ -297,7 +297,7 @@ end
 
 #     chain = Chain(kernel, μ)
 
-#     dyniterate(iter::TimeLift, ::Nothing) = dyniterate(iter, 0=>nothing) 
+#     dyniterate(iter::TimeLift, ::Nothing) = dyniterate(iter, 0=>nothing)
 #     tr1 = trace(TimeLift(chain), nothing, u -> u[1] > 15)
 #     tr2 = trace(TimeLift(rand(Random.GLOBAL_RNG, chain)), nothing, u -> u[1] > 15)
 #     collect(Iterators.take(chain, 10))
@@ -348,8 +348,8 @@ end
     # NOTE: The `test_broken` below are mostly because of the change to `Affine`.
     # For example, `Normal{(:μ,:σ)}` is now `Affine{(:μ,:σ), Normal{()}}`.
     # The problem is not really with these measures, but with the tests
-    # themselves. 
-    # 
+    # themselves.
+    #
     # We should instead probably be doing e.g.
     # `D = typeof(Normal(μ=0.3, σ=4.1))`
 
@@ -652,11 +652,29 @@ end
     end
 
     x = rand(d)
-    
+
     @test logdensityof(d, x) isa Real
 end
 
 @testset "Distributions.jl cdf" begin
     @test cdf(Normal(0, 1), 0) == 0.5
     @test cdf.((Normal(0, 1),), [0, 0]) == [0.5, 0.5]
 end
+
+@testset "pairwise normal-normal" begin
+    n0 = Normal(0, 0)
+    n1 = Normal(0.0, 1.0)
+    n2 = Normal(1.0, 2.0)
+    @test MeasureTheory.convolve(n1, n0) == n1
+    @test MeasureTheory.convolve(n1, n2) == MeasureTheory.convolve(n2, n1)
+
+    standard_normal_kernel = MeasureTheory.kernel(Normal{(:μ,)}, μ=identity)
+    σs = [0.0, 1.0, 2.0]
+    normal_kernels = [MeasureTheory.kernel(Normal{(:μ,:σ)}, μ=identity, σ=σ) for σ in σs]
+
+    @test (n1 ↣ standard_normal_kernel) == Normal(0.0, sqrt(2))
+    @test (n0 ↣ standard_normal_kernel) == n1
+
+    @test (n2 ↣ normal_kernels[1]) == n2
+    @test (n2 ↣ normal_kernels[2]) == (n2 ↣ standard_normal_kernel)
+end