Repository: biaslab/ReactiveMP.jl (julia)

Issue:
Outgoing BP rules for :var interface of NormalMeanVariance node.
### Discussed in https://github.com/biaslab/ReactiveMP.jl/discussions/95

Hi guys,
(@ismailsenoz @semihakbayrak @bvdmitri)
How sure we are about the BP message for `:var` interface that is associated with variance?

We had a look a few days ago with @bartvanerp and it seemed to us that the normalization of the message is incorrect.
I see there are tests for these rules, but I am still concerned about some missing factors in:
```julia
ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_μ_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
```


There is an issue with BP message from `GaussianMeanVariance` node. The correct version when both `:mean` and `:out` inbounds messages are Gaussians should be:
```
(x) -> -log(m_μ_cov + m_out_cov + x) / 2 - 1 / (2*(m_μ_cov + m_out_cov + x)) * (m_out_mean - m_μ_mean) ^ 2)
```

Pull request description:
Fixes #96. @albertpod, @bartvanerp, @ismailsenoz please double check.

Patch:
diff --git a/src/rules/normal_mean_variance/var.jl b/src/rules/normal_mean_variance/var.jl
index cac0e892b..8ac3b2a95 100644
--- a/src/rules/normal_mean_variance/var.jl
+++ b/src/rules/normal_mean_variance/var.jl
@@ -5,19 +5,19 @@
     m_out_mean        = mean(m_out)
     m_μ_mean, m_μ_cov = mean_cov(m_μ)
     
-    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_μ_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
+    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_μ_cov + x) / 2 - 1 / (2*(x + m_μ_cov)) * (m_out_mean - m_μ_mean) ^ 2)
 end
 
 @rule NormalMeanVariance(:v, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_μ::PointMass) = begin 
     m_out_mean, m_out_cov = mean_cov(m_out)
     m_μ_mean              = mean(m_μ)
 
-    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
+    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + x) / 2 - 1 / (2*(x + m_out_cov)) * (m_out_mean - m_μ_mean) ^ 2)
 end
 
 @rule NormalMeanVariance(:v, Marginalisation) (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily) = begin 
     m_out_mean, m_out_cov = mean_cov(m_out)
     m_μ_mean, m_μ_cov     = mean_cov(m_μ)
 
-    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + m_μ_cov + x) / 2 - 1 / (2x) * (m_out_mean - m_μ_mean) ^ 2)
+    return ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(m_out_cov + m_μ_cov + x) / 2 - 1 / (2*(m_out_cov + m_μ_cov + x)) * (m_out_mean - m_μ_mean) ^ 2)
 end


Test patch:
diff --git a/test/rules/normal_mean_variance/test_var.jl b/test/rules/normal_mean_variance/test_var.jl
index 4caff03f3..9ef805da8 100644
--- a/test/rules/normal_mean_variance/test_var.jl
+++ b/test/rules/normal_mean_variance/test_var.jl
@@ -12,10 +12,10 @@ import ReactiveMP: @test_rules
     @testset "Belief Propagation: (m_out::PointMass, m_μ::UnivariateNormalDistributionsFamily)" begin
 
         @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
-            (input = (m_out = PointMass(2.0), m_μ = NormalMeanVariance(0.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / 2x)),
-            (input = (m_out = PointMass(0.5), m_μ = NormalMeanVariance(1.0, 0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / 2x)),
-            (input = (m_out = PointMass(-3.5), m_μ = NormalMeanVariance(0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
-            (input = (m_out = PointMass(3.5), m_μ = NormalMeanVariance(-0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
+            (input = (m_out = PointMass(2.0), m_μ = NormalMeanVariance(0.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / (2*(1.0 + x)))),
+            (input = (m_out = PointMass(0.5), m_μ = NormalMeanVariance(1.0, 0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / (2*(0.5 + x)))),
+            (input = (m_out = PointMass(-3.5), m_μ = NormalMeanVariance(0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / (2*(2.0 + x)))),
+            (input = (m_out = PointMass(3.5), m_μ = NormalMeanVariance(-0.5, 2.0)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / (2*(2.0 + x)))),
         ]
 
     end
@@ -23,10 +23,10 @@ import ReactiveMP: @test_rules
     @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_μ::PointMass)" begin
 
         @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
-            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = PointMass(2.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / 2x)),
-            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = PointMass(0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / 2x)),
-            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = PointMass(-3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
-            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = PointMass(3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / 2x)),
+            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = PointMass(2.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(1.0 + x) / 2 - 4 / (2*(1.0 + x)))),
+            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = PointMass(0.5)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(0.5 + x) / 2 - 0.25 / (2*(0.5 + x)))),
+            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = PointMass(-3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / (2*(2.0 + x)))),
+            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = PointMass(3.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 16 / (2*(2.0 + x)))),
         ]
 
     end
@@ -34,10 +34,10 @@ import ReactiveMP: @test_rules
     @testset "Belief Propagation: (m_out::UnivariateNormalDistributionsFamily, m_μ::UnivariateNormalDistributionsFamily)" begin
 
         @test_rules [ with_float_conversions = true ] NormalMeanVariance(:v, Marginalisation) [
-            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = NormalMeanVariance(-2.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 4 / 2x)),
-            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = NormalMeanVariance(0.0, 2.0)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
-            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = NormalMeanVariance(1.5, 0.5)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
-            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = NormalMeanVariance(-1.5, 0.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / 2x)),
+            (input = (m_out = NormalMeanVariance(0.0, 1.0), m_μ = NormalMeanVariance(-2.0, 1.0)),  output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.0 + x) / 2 - 4 / (2*(2.0 + x)))),
+            (input = (m_out = NormalMeanVariance(1.0, 0.5), m_μ = NormalMeanVariance(0.0, 2.0)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / (2*(2.5 + x)))),
+            (input = (m_out = NormalMeanVariance(0.5, 2.0), m_μ = NormalMeanVariance(1.5, 0.5)),   output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / (2*(2.5 + x)))),
+            (input = (m_out = NormalMeanVariance(-0.5, 2.0), m_μ = NormalMeanVariance(-1.5, 0.5)), output = ContinuousUnivariateLogPdf(DomainSets.HalfLine(), (x) -> -log(2.5 + x) / 2 - 1 / (2*(2.5 + x)))),
         ]
 
     end