Merge pull request #418 from ReactiveBayes/multinomial_polya_model

Multinomial polya model

Project.toml

@@ -1,7 +1,7 @@
 name = "RxInfer"
 uuid = "86711068-29c9-4ff7-b620-ae75d7495b3d"
 authors = ["Bagaev Dmitry <d.v.bagaev@tue.nl> and contributors"]
-version = "4.0.1"
+version = "4.1.0"
 
 [deps]
 BayesBase = "b4ee3484-f114-42fe-b91c-797d54a0c67e"
@@ -56,7 +56,7 @@ Preferences = "1.4.3"
 PrettyTables = "2"
 ProgressMeter = "1.0.0"
 Random = "1.9"
-ReactiveMP = "~5.0.0"
+ReactiveMP = "~5.1.0"
 Reexport = "1.2.0"
 Rocket = "1.8.0"
 Static = "0.8.10, 1"

codemeta.json

@@ -9,12 +9,12 @@
     "downloadUrl": "/~https://github.com/reactivebayes/RxInfer.jl/releases",
     "issueTracker": "/~https://github.com/reactivebayes/RxInfer.jl/issues",
     "name": "RxInfer.jl",
-    "version": "4.0.1",
+    "version": "4.1.0",
     "description": "Julia package for automated, scalable and efficient Bayesian inference on factor graphs with reactive message passing. ",
     "applicationCategory": "Statistics",
     "developmentStatus": "active",
     "readme": "https://rxinfer.ml",
-    "softwareVersion": "4.0.1",
+    "softwareVersion": "4.1.0",
     "keywords": [
         "Bayesian inference",
         "message passing",

test/models/regression/multinomialreg_tests.jl

@@ -0,0 +1,81 @@
+@testitem "Multinomial regression with MultinomialPolya (offline inference) node" begin
+    using BenchmarkTools, Plots, Distributions, LinearAlgebra, StableRNGs, ExponentialFamily.LogExpFunctions
+
+    include(joinpath(@__DIR__, "..", "..", "utiltests.jl"))
+
+    N = 20
+    k = 10
+    nsamples = 1000
+    X, ψ, p = generate_multinomial_data(; N = N, k = k, nsamples = nsamples)
+
+    @model function multinomial_model(y, N, ξ_ψ, W_ψ)
+        ψ ~ MvNormalWeightedMeanPrecision(ξ_ψ, W_ψ)
+        for i in eachindex(y)
+            y[i] ~ MultinomialPolya(N, ψ) where {dependencies = RequireMessageFunctionalDependencies(ψ = MvNormalWeightedMeanPrecision(ξ_ψ, W_ψ))}
+        end
+    end
+
+    result = infer(
+        model = multinomial_model(ξ_ψ = zeros(k - 1), W_ψ = rand(Wishart(k, diageye(k - 1))), N = N),
+        data = (y = X,),
+        iterations = 100,
+        free_energy = true,
+        showprogress = false,
+        returnvars = KeepLast(),
+        options = (limit_stack_depth = 100,)
+    )
+
+    m = mean(result.posteriors[:ψ])
+    pest = logistic_stic_breaking(m)
+
+    mse = mean((pest - p) .^ 2)
+    @test mse < 2e-5
+
+    @test result.free_energy[end] < result.free_energy[1]
+    @test result.free_energy[end] <= result.free_energy[end - 1]
+    @test abs(result.free_energy[end - 1] - result.free_energy[end]) < 1e-8
+end
+
+@testitem "Multinomial regression - online inference" begin
+    using BenchmarkTools, Plots, Distributions, LinearAlgebra, StableRNGs, ExponentialFamily.LogExpFunctions
+
+    include(joinpath(@__DIR__, "..", "..", "utiltests.jl"))
+
+    N = 50
+    k = 40
+    nsamples = 5000
+    X, ψ, p = generate_multinomial_data(; N = N, k = k, nsamples = nsamples)
+
+    @model function multinomial_model(y, N, ξ_ψ, W_ψ, k)
+        ψ ~ MvNormalWeightedMeanPrecision(ξ_ψ, W_ψ)
+        y ~ MultinomialPolya(N, ψ) where {dependencies = RequireMessageFunctionalDependencies(ψ = MvNormalWeightedMeanPrecision(zeros(k - 1), diageye(k - 1)))}
+    end
+
+    @autoupdates function auto()
+        ξ_ψ, W_ψ = weightedmean_precision(q(ψ))
+    end
+    init = @initialization begin
+        q(ψ) = MvNormalWeightedMeanPrecision(zeros(k - 1), rand(Wishart(k, diageye(k - 1))))
+    end
+
+    result = infer(
+        model = multinomial_model(N = N, k = k),
+        data = (y = X,),
+        initialization = init,
+        iterations = 1,
+        autoupdates = auto(),
+        keephistory = length(X),
+        free_energy = true,
+        showprogress = false
+    )
+
+    m = result.history[:ψ][end]
+
+    pest = logistic_stic_breaking(mean(m))
+    mse = mean((pest - p) .^ 2)
+    @test mse < 1e-3
+
+    @test result.free_energy_final_only_history[end] < result.free_energy_final_only_history[1]
+    #Free energy over time decreases in a noisy way. It is not a monotonic decrease. 
+
+end

test/utiltests.jl

@@ -79,3 +79,24 @@ macro test_expression_generating(lhs, rhs)
         end
     )
 end
+
+function generate_multinomial_data(rng = StableRNG(123); N = 3, k = 3, nsamples = 5000)
+    ψ = randn(rng, k)
+    p = ReactiveMP.softmax(ψ)
+
+    X = rand(rng, Multinomial(N, p), nsamples)
+    X = [X[:, i] for i in axes(X, 2)]
+    return X, ψ, p
+end
+
+function logistic_stic_breaking(m)
+    Km1 = length(m)
+
+    p = Array{Float64}(undef, Km1 + 1)
+    p[1] = logistic(m[1])
+    for i in 2:Km1
+        p[i] = logistic(m[i]) * (1 - sum(p[1:(i - 1)]))
+    end
+    p[end] = 1 - sum(p[1:(end - 1)])
+    return p
+end