From 8b06dee80beda9b239feb26eca8dc29618ff8ff9 Mon Sep 17 00:00:00 2001 From: Xianda Sun Date: Mon, 7 Oct 2024 05:58:07 +0100 Subject: [PATCH 1/5] simplify log density tests --- test/logp_tests/test_logp.jl | 294 +++++++++++++++++++++++++++++++++++ 1 file changed, 294 insertions(+) diff --git a/test/logp_tests/test_logp.jl b/test/logp_tests/test_logp.jl index a6a34aec8..931bd49bd 100644 --- a/test/logp_tests/test_logp.jl +++ b/test/logp_tests/test_logp.jl @@ -14,3 +14,297 @@ using DynamicPPL: DynamicPPL, getlogp, settrans!!, SimpleVarInfo, @model end end end + +function _logjoint(model::JuliaBUGS.BUGSModel) + return JuliaBUGS.evaluate!!(model, JuliaBUGS.DefaultContext())[2] +end + +@testset "Log density" begin + @testset "Log density of distributions" begin + @testset "dbin (Binomial)" begin + dist = dbin(0.1, 10) + b = Bijectors.bijector(dist) + test_θ_transformed = 10 + test_θ = Bijectors.inverse(b)(test_θ_transformed) + + model_def = @bugs begin + a ~ dbin(0.1, 10) + end + transformed_model = compile(model_def, NamedTuple(), (a=test_θ,)) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + + reference_logp_untransformed = logpdf(dist, test_θ) + reference_logp_transformed = + logpdf(dist, test_θ) + + logabsdetjac(Bijectors.inverse(b), test_θ_transformed) + + # the bijector of dbin is the identity, so the log density should be the same + @test _logjoint(untransformed_model) ≈ reference_logp_untransformed rtol = 1E-6 + @test _logjoint(transformed_model) ≈ reference_logp_transformed rtol = 1E-6 + end + + @testset "dgamma (Gamma)" begin + dist = dgamma(0.001, 0.001) + b = Bijectors.bijector(dist) + test_θ_transformed = 10 + test_θ = Bijectors.inverse(b)(test_θ_transformed) + + model_def = @bugs begin + a ~ dgamma(0.001, 0.001) + end + transformed_model = compile(model_def, NamedTuple(), (a=test_θ,)) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + + reference_logp_untransformed = logpdf(dist, test_θ) + reference_logp_transformed = + logpdf(dist, test_θ) + + logabsdetjac(Bijectors.inverse(b), test_θ_transformed) + + @test _logjoint(untransformed_model) ≈ reference_logp_untransformed rtol = 1E-6 + @test _logjoint(transformed_model) ≈ reference_logp_transformed rtol = 1E-6 + end + + @testset "ddirich (Dirichlet)" begin + # create valid test input + alpha = rand(10) + dist = ddirich(alpha) + b = Bijectors.bijector(dist) + test_θ_transformed = rand(9) + test_θ = Bijectors.inverse(b)(test_θ_transformed) + + reference_logp_untransformed = logpdf(dist, test_θ) + reference_logp_transformed = + logpdf(dist, test_θ) + + logabsdetjac(Bijectors.inverse(b), test_θ_transformed) + + model_def = @bugs begin + x[1:10] ~ ddirich(alpha[1:10]) + end + transformed_model = compile(model_def, (alpha=alpha,), (x=test_θ,)) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + + @test _logjoint(untransformed_model) ≈ reference_logp_untransformed rtol = 1E-6 + @test _logjoint(transformed_model) ≈ reference_logp_transformed rtol = 1E-6 + end + + @testset "dwish (Wishart)" begin + # create valid test input + scale_matrix = randn(10, 10) + scale_matrix = scale_matrix * transpose(scale_matrix) # Ensuring positive-definiteness + degrees_of_freedom = 12 + + dist = dwish(scale_matrix, degrees_of_freedom) + b = Bijectors.bijector(dist) + test_θ_transformed = rand(55) + test_θ = Bijectors.inverse(b)(test_θ_transformed) + + reference_logp_untransformed = logpdf(dist, test_θ) + reference_logp_transformed = + logpdf(dist, test_θ) + + logabsdetjac(Bijectors.inverse(b), test_θ_transformed) + + model_def = @bugs begin + x[1:10, 1:10] ~ dwish(scale_matrix[:, :], degrees_of_freedom) + end + transformed_model = compile( + model_def, + (degrees_of_freedom=degrees_of_freedom, scale_matrix=scale_matrix), + (x=test_θ,), + ) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + + @test _logjoint(untransformed_model) ≈ reference_logp_untransformed rtol = 1E-6 + @test _logjoint(transformed_model) ≈ reference_logp_transformed rtol = 1E-6 + end + + @testset "lkj (LKJ)" begin + dist = LKJ(10, 0.5) + b = Bijectors.bijector(dist) + test_θ_transformed = rand(45) + test_θ = Bijectors.inverse(b)(test_θ_transformed) + + reference_logp_untransformed = logpdf(dist, test_θ) + reference_logp_transformed = + logpdf(dist, test_θ) + + logabsdetjac(Bijectors.inverse(b), test_θ_transformed) + + model_def = @bugs begin + x[1:10, 1:10] ~ LKJ(10, 0.5) + end + transformed_model = compile(model_def, NamedTuple(), (x=test_θ,)) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + + @test LogDensityProblems.dimension(untransformed_model) == 100 + @test LogDensityProblems.dimension(transformed_model) == 45 + + @test _logjoint(untransformed_model) ≈ reference_logp_untransformed rtol = 1E-6 + @test _logjoint(transformed_model) ≈ reference_logp_transformed rtol = 1E-6 + end + end +end + +@testset "Log density of BUGS models" begin + @testset "rats" begin + (; model_def, data, inits) = JuliaBUGS.BUGSExamples.VOLUME_1.rats + transformed_model = compile(model_def, data, inits) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + @test _logjoint(untransformed_model) ≈ -174029.38703951868 rtol = 1E-6 + @test _logjoint(transformed_model) ≈ -174029.38703951868 rtol = 1E-6 + end + + @testset "blockers" begin + (; model_def, data, inits) = JuliaBUGS.BUGSExamples.VOLUME_1.blockers + transformed_model = compile(model_def, data, inits) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + @test _logjoint(untransformed_model) ≈ -8418.416388326123 rtol = 1E-6 + @test _logjoint(transformed_model) ≈ -8418.416388326123 rtol = 1E-6 + end + + @testset "bones" begin + (; model_def, data, inits) = JuliaBUGS.BUGSExamples.VOLUME_1.bones + transformed_model = compile(model_def, data, inits) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + @test _logjoint(untransformed_model) ≈ -161.6492002285034 rtol = 1E-6 + @test _logjoint(transformed_model) ≈ -161.6492002285034 rtol = 1E-6 + end + + @testset "dogs" begin + (; model_def, data, inits) = JuliaBUGS.BUGSExamples.VOLUME_1.dogs + transformed_model = compile(model_def, data, inits) + untransformed_model = JuliaBUGS.settrans(transformed_model, false) + @test _logjoint(untransformed_model) ≈ -1243.188922285352 rtol = 1E-6 + @test _logjoint(transformed_model) ≈ -1243.3996613167667 rtol = 1E-6 + end +end + +## transcribed BUGS models in DynamicPPL + +# # rats +# @model function rats(Y, x, xbar, N, T) +# var"tau.c" ~ dgamma(0.001, 0.001) +# sigma = 1 / sqrt(var"tau.c") + +# var"alpha.c" ~ dnorm(0.0, 1.0E-6) +# var"alpha.tau" ~ dgamma(0.001, 0.001) + +# var"beta.c" ~ dnorm(0.0, 1.0E-6) +# var"beta.tau" ~ dgamma(0.001, 0.001) + +# alpha0 = var"alpha.c" - xbar * var"beta.c" + +# alpha = Vector{Real}(undef, N) +# beta = Vector{Real}(undef, N) + +# for i in 1:N +# alpha[i] ~ dnorm(var"alpha.c", var"alpha.tau") +# beta[i] ~ dnorm(var"beta.c", var"beta.tau") + +# for j in 1:T +# mu = alpha[i] + beta[i] * (x[j] - xbar) +# Y[i, j] ~ dnorm(mu, var"tau.c") +# end +# end + +# return sigma, alpha0 +# end + +# (; N, T, x, xbar, Y) = data +# model = rats(Y, x, xbar, N, T) + +# # blockers +# @model function blockers(rc, rt, nc, nt, Num) +# d ~ dnorm(0.0, 1.0E-6) +# tau ~ dgamma(0.001, 0.001) + +# mu = Vector{Real}(undef, Num) +# delta = Vector{Real}(undef, Num) +# pc = Vector{Real}(undef, Num) +# pt = Vector{Real}(undef, Num) + +# for i in 1:Num +# mu[i] ~ dnorm(0.0, 1.0E-5) +# delta[i] ~ dnorm(d, tau) + +# pc[i] = logistic(mu[i]) +# pt[i] = logistic(mu[i] + delta[i]) + +# rc[i] ~ dbin(pc[i], nc[i]) +# rt[i] ~ dbin(pt[i], nt[i]) +# end + +# var"delta.new" ~ dnorm(d, tau) +# sigma = 1 / sqrt(tau) + +# return sigma +# end + +# (; rt, nt, rc, nc, Num) = data +# model = blockers(rc, rt, nc, nt, Num) + +# # bones +# @model function bones(grade, nChild, nInd, ncat, gamma, delta) +# theta = Vector{Real}(undef, nChild) +# Q = Array{Real}(undef, nChild, nInd, maximum(ncat)) +# p = Array{Real}(undef, nChild, nInd, maximum(ncat)) +# cumulative_grade = Array{Real}(undef, nChild, nInd) + +# for i in 1:nChild +# theta[i] ~ dnorm(0.0, 0.001) + +# for j in 1:nInd +# for k in 1:(ncat[j] - 1) +# Q[i, j, k] = logistic(delta[j] * (theta[i] - gamma[j, k])) +# end +# end + +# for j in 1:nInd +# p[i, j, 1] = 1 - Q[i, j, 1] + +# for k in 2:(ncat[j] - 1) +# p[i, j, k] = Q[i, j, k - 1] - Q[i, j, k] +# end + +# p[i, j, ncat[j]] = Q[i, j, ncat[j] - 1] +# grade[i, j] ~ dcat(p[i, j, 1:ncat[j]]) +# end +# end +# end + +# (; grade, nChild, nInd, ncat, gamma, delta) = data +# model = bones(grade, nChild, nInd, ncat, gamma, delta) + +# # dogs + +# @model function dogs(Dogs, Trials, Y, y) +# # Initialize matrices +# xa = zeros(Dogs, Trials) +# xs = zeros(Dogs, Trials) +# p = zeros(Dogs, Trials) + +# # Flat priors for alpha and beta, restricted to (-∞, -0.00001) +# alpha ~ dunif(-10, -1.0e-5) +# beta ~ dunif(-10, -1.0e-5) + +# for i in 1:Dogs +# xa[i, 1] = 0 +# xs[i, 1] = 0 +# p[i, 1] = 0 + +# for j in 2:Trials +# xa[i, j] = sum(Y[i, 1:(j - 1)]) +# xs[i, j] = j - 1 - xa[i, j] +# p[i, j] = exp(alpha * xa[i, j] + beta * xs[i, j]) +# # The Bernoulli likelihood +# y[i, j] ~ dbern(p[i, j]) +# end +# end + +# # Transformation to positive values +# A = exp(alpha) +# B = exp(beta) + +# return A, B +# end + +# (; Dogs, Trials, Y) = data +# model = dogs(Dogs, Trials, Y, 1 .- Y) From 2007dda4ce0b3b951df3e6f1d9a041aafc134bb1 Mon Sep 17 00:00:00 2001 From: Xianda Sun Date: Mon, 7 Oct 2024 05:59:02 +0100 Subject: [PATCH 2/5] move logp test file --- test/runtests.jl | 2 +- test/{logp_tests => }/test_logp.jl | 0 2 files changed, 1 insertion(+), 1 deletion(-) rename test/{logp_tests => }/test_logp.jl (100%) diff --git a/test/runtests.jl b/test/runtests.jl index 37cb50f4a..b62ce06b0 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ if test_group == "compilation" || test_group == "all" include("bugs_primitives.jl") include("compile.jl") end - include("logp_tests/test_logp.jl") + include("test_logp.jl") end if test_group == "gibbs" || test_group == "all" diff --git a/test/logp_tests/test_logp.jl b/test/test_logp.jl similarity index 100% rename from test/logp_tests/test_logp.jl rename to test/test_logp.jl From 87c36e74fdc3179bc655ad0952354f4aa96603ea Mon Sep 17 00:00:00 2001 From: Xianda Sun Date: Mon, 7 Oct 2024 06:01:15 +0100 Subject: [PATCH 3/5] rename file --- test/{test_logp.jl => log_density.jl} | 14 +---- test/logp_tests/BUGS_models/blockers.jl | 58 ----------------- test/logp_tests/BUGS_models/bones.jl | 60 ------------------ test/logp_tests/BUGS_models/dogs.jl | 56 ----------------- test/logp_tests/BUGS_models/rats.jl | 63 ------------------- .../single_distribution_models/binomial.jl | 35 ----------- .../single_distribution_models/ddirich.jl | 46 -------------- .../single_distribution_models/dwish.jl | 53 ---------------- .../single_distribution_models/gamma.jl | 36 ----------- .../single_distribution_models/lkj.jl | 44 ------------- test/runtests.jl | 2 +- 11 files changed, 2 insertions(+), 465 deletions(-) rename test/{test_logp.jl => log_density.jl} (96%) delete mode 100644 test/logp_tests/BUGS_models/blockers.jl delete mode 100644 test/logp_tests/BUGS_models/bones.jl delete mode 100644 test/logp_tests/BUGS_models/dogs.jl delete mode 100644 test/logp_tests/BUGS_models/rats.jl delete mode 100644 test/logp_tests/single_distribution_models/binomial.jl delete mode 100644 test/logp_tests/single_distribution_models/ddirich.jl delete mode 100644 test/logp_tests/single_distribution_models/dwish.jl delete mode 100644 test/logp_tests/single_distribution_models/gamma.jl delete mode 100644 test/logp_tests/single_distribution_models/lkj.jl diff --git a/test/test_logp.jl b/test/log_density.jl similarity index 96% rename from test/test_logp.jl rename to test/log_density.jl index 931bd49bd..3f1ed3489 100644 --- a/test/test_logp.jl +++ b/test/log_density.jl @@ -2,19 +2,7 @@ using JuliaBUGS: BUGSGraph, DefaultContext, evaluate!!, get_params_varinfo, LogDensityContext using DynamicPPL: DynamicPPL, getlogp, settrans!!, SimpleVarInfo, @model -@testset "Log joint probability" begin - @testset "Single distribution models" begin - @testset "$s" for s in [:binomial, :gamma, :lkj, :dwish, :ddirich] - include("single_distribution_models/$s.jl") - end - end - @testset "BUGS models" begin - @testset "$s" for s in [:blockers, :bones, :dogs, :rats] - include("BUGS_models/$s.jl") - end - end -end - +# TODO: make this available in JuliaBUGS function _logjoint(model::JuliaBUGS.BUGSModel) return JuliaBUGS.evaluate!!(model, JuliaBUGS.DefaultContext())[2] end diff --git a/test/logp_tests/BUGS_models/blockers.jl b/test/logp_tests/BUGS_models/blockers.jl deleted file mode 100644 index eda42f277..000000000 --- a/test/logp_tests/BUGS_models/blockers.jl +++ /dev/null @@ -1,58 +0,0 @@ -bugs_model_def = JuliaBUGS.BUGSExamples.VOLUME_1.blockers.model_def -data = JuliaBUGS.BUGSExamples.VOLUME_1.blockers.data -inits = JuliaBUGS.BUGSExamples.VOLUME_1.blockers.inits - -bugs_model = compile(bugs_model_def, data, inits) -vi = bugs_model.varinfo - -@model function blockers(rc, rt, nc, nt, Num) - d ~ dnorm(0.0, 1.0E-6) - tau ~ dgamma(0.001, 0.001) - - mu = Vector{Real}(undef, Num) - delta = Vector{Real}(undef, Num) - pc = Vector{Real}(undef, Num) - pt = Vector{Real}(undef, Num) - - for i in 1:Num - mu[i] ~ dnorm(0.0, 1.0E-5) - delta[i] ~ dnorm(d, tau) - - pc[i] = logistic(mu[i]) - pt[i] = logistic(mu[i] + delta[i]) - - rc[i] ~ dbin(pc[i], nc[i]) - rt[i] ~ dbin(pt[i], nt[i]) - end - - var"delta.new" ~ dnorm(d, tau) - sigma = 1 / sqrt(tau) - - return sigma -end - -(; rt, nt, rc, nc, Num) = data -dppl_model = blockers(rc, rt, nc, nt, Num) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -params_vi = JuliaBUGS.get_params_varinfo(bugs_model, vi) -# test if JuliaBUGS and DynamicPPL agree on parameters in the model -@test keys( - DynamicPPL.evaluate!!( - dppl_model, SimpleVarInfo(Dict{VarName,Any}()), DynamicPPL.SamplingContext() - )[2], -) == keys(params_vi) - -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, DynamicPPL.settrans!!(vi, false), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ -8418.416388 rtol = 1E-6 -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, get_params_varinfo(bugs_model), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/logp_tests/BUGS_models/bones.jl b/test/logp_tests/BUGS_models/bones.jl deleted file mode 100644 index 40a709901..000000000 --- a/test/logp_tests/BUGS_models/bones.jl +++ /dev/null @@ -1,60 +0,0 @@ -bugs_model_def = JuliaBUGS.BUGSExamples.VOLUME_1[:bones].model_def -data = JuliaBUGS.BUGSExamples.VOLUME_1[:bones].data -inits = JuliaBUGS.BUGSExamples.VOLUME_1[:bones].inits - -bugs_model = compile(bugs_model_def, data, inits) -vi = bugs_model.varinfo - -@model function bones(grade, nChild, nInd, ncat, gamma, delta) - theta = Vector{Real}(undef, nChild) - Q = Array{Real}(undef, nChild, nInd, maximum(ncat)) - p = Array{Real}(undef, nChild, nInd, maximum(ncat)) - cumulative_grade = Array{Real}(undef, nChild, nInd) - - for i in 1:nChild - theta[i] ~ dnorm(0.0, 0.001) - - for j in 1:nInd - for k in 1:(ncat[j] - 1) - Q[i, j, k] = logistic(delta[j] * (theta[i] - gamma[j, k])) - end - end - - for j in 1:nInd - p[i, j, 1] = 1 - Q[i, j, 1] - - for k in 2:(ncat[j] - 1) - p[i, j, k] = Q[i, j, k - 1] - Q[i, j, k] - end - - p[i, j, ncat[j]] = Q[i, j, ncat[j] - 1] - grade[i, j] ~ dcat(p[i, j, 1:ncat[j]]) - end - end -end - -(; grade, nChild, nInd, ncat, gamma, delta) = data -dppl_model = bones(grade, nChild, nInd, ncat, gamma, delta) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -params_vi = JuliaBUGS.get_params_varinfo(bugs_model, vi) -# test if JuliaBUGS and DynamicPPL agree on parameters in the model -@test keys( - DynamicPPL.evaluate!!( - dppl_model, SimpleVarInfo(Dict{VarName,Any}()), DynamicPPL.SamplingContext() - )[2], -) == keys(params_vi) - -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, DynamicPPL.settrans!!(vi, false), DynamicPPL.DefaultContext() - )[2].logp -# ! ProbPALA compile error -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, get_params_varinfo(bugs_model), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/logp_tests/BUGS_models/dogs.jl b/test/logp_tests/BUGS_models/dogs.jl deleted file mode 100644 index e60b812c2..000000000 --- a/test/logp_tests/BUGS_models/dogs.jl +++ /dev/null @@ -1,56 +0,0 @@ -bugs_model_def = JuliaBUGS.BUGSExamples.VOLUME_1[:dogs].model_def -data = JuliaBUGS.BUGSExamples.VOLUME_1[:dogs].data -inits = JuliaBUGS.BUGSExamples.VOLUME_1[:dogs].inits - -bugs_model = compile(bugs_model_def, data, inits) -vi = bugs_model.varinfo - -@model function dogs(Dogs, Trials, Y, y) - # Initialize matrices - xa = zeros(Dogs, Trials) - xs = zeros(Dogs, Trials) - p = zeros(Dogs, Trials) - - # Flat priors for alpha and beta, restricted to (-∞, -0.00001) - alpha ~ dunif(-10, -1.0e-5) - beta ~ dunif(-10, -1.0e-5) - - for i in 1:Dogs - xa[i, 1] = 0 - xs[i, 1] = 0 - p[i, 1] = 0 - - for j in 2:Trials - xa[i, j] = sum(Y[i, 1:(j - 1)]) - xs[i, j] = j - 1 - xa[i, j] - p[i, j] = exp(alpha * xa[i, j] + beta * xs[i, j]) - # The Bernoulli likelihood - y[i, j] ~ dbern(p[i, j]) - end - end - - # Transformation to positive values - A = exp(alpha) - B = exp(beta) - - return A, B -end - -(; Dogs, Trials, Y) = data -dppl_model = dogs(Dogs, Trials, Y, 1 .- Y) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] - -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, DynamicPPL.settrans!!(vi, false), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ -1243.188922 rtol = 1E-6 # reference value from ProbPALA -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, get_params_varinfo(bugs_model), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/logp_tests/BUGS_models/rats.jl b/test/logp_tests/BUGS_models/rats.jl deleted file mode 100644 index 716a53035..000000000 --- a/test/logp_tests/BUGS_models/rats.jl +++ /dev/null @@ -1,63 +0,0 @@ -model_def = JuliaBUGS.BUGSExamples.VOLUME_1.rats.model_def -data = JuliaBUGS.BUGSExamples.VOLUME_1.rats.data -inits = JuliaBUGS.BUGSExamples.VOLUME_1.rats.inits - -bugs_model = compile(model_def, data, inits); -vi = bugs_model.varinfo - -@model function rats(Y, x, xbar, N, T) - var"tau.c" ~ dgamma(0.001, 0.001) - sigma = 1 / sqrt(var"tau.c") - - var"alpha.c" ~ dnorm(0.0, 1.0E-6) - var"alpha.tau" ~ dgamma(0.001, 0.001) - - var"beta.c" ~ dnorm(0.0, 1.0E-6) - var"beta.tau" ~ dgamma(0.001, 0.001) - - alpha0 = var"alpha.c" - xbar * var"beta.c" - - alpha = Vector{Real}(undef, N) - beta = Vector{Real}(undef, N) - - for i in 1:N - alpha[i] ~ dnorm(var"alpha.c", var"alpha.tau") - beta[i] ~ dnorm(var"beta.c", var"beta.tau") - - for j in 1:T - mu = alpha[i] + beta[i] * (x[j] - xbar) - Y[i, j] ~ dnorm(mu, var"tau.c") - end - end - - return sigma, alpha0 -end -(; N, T, x, xbar, Y) = data -dppl_model = rats(Y, x, xbar, N, T) - -bugs_model = JuliaBUGS.settrans(bugs_model, false) -bugs_logp = JuliaBUGS.evaluate!!(bugs_model, DefaultContext())[2] -params_vi = JuliaBUGS.get_params_varinfo(bugs_model, vi) -# test if JuliaBUGS and DynamicPPL agree on parameters in the model -@test keys( - DynamicPPL.evaluate!!( - dppl_model, SimpleVarInfo(Dict{VarName,Any}()), DynamicPPL.SamplingContext() - )[2], -) == keys(params_vi) - -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, DynamicPPL.settrans!!(vi, false), DynamicPPL.DefaultContext() - )[2].logp -@test bugs_logp ≈ -174029.387 rtol = 1E-6 # reference value from ProbPALA -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -dppl_logp = - DynamicPPL.evaluate!!( - dppl_model, get_params_varinfo(bugs_model), DynamicPPL.DefaultContext() - )[2].logp -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -@test bugs_model.untransformed_param_length == - LogDensityProblems.dimension(DynamicPPL.LogDensityFunction(dppl_model)) diff --git a/test/logp_tests/single_distribution_models/binomial.jl b/test/logp_tests/single_distribution_models/binomial.jl deleted file mode 100644 index e0620047f..000000000 --- a/test/logp_tests/single_distribution_models/binomial.jl +++ /dev/null @@ -1,35 +0,0 @@ -model_def = @bugs begin - a ~ dbin(0.1, 10) -end - -bugs_model = compile(model_def, NamedTuple(), (a=10,)) -vi = bugs_model.varinfo - -@model function dppl_gamma_model() - return a ~ dbin(0.1, 10) -end - -dppl_model = dppl_gamma_model() - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -params_vi = JuliaBUGS.get_params_varinfo(bugs_model, vi) -# test if JuliaBUGS and DynamicPPL agree on parameters in the model -@test keys( - DynamicPPL.evaluate!!( - dppl_model, SimpleVarInfo(Dict{VarName,Any}()), DynamicPPL.SamplingContext() - )[2], -) == keys(params_vi) - -p = DynamicPPL.LogDensityFunction(dppl_model) -t_p = DynamicPPL.LogDensityFunction( - dppl_model, - DynamicPPL.link!!(SimpleVarInfo(dppl_model), dppl_model), - DynamicPPL.DefaultContext(), -) - -@test bugs_logp ≈ LogDensityProblems.logdensity(p, [10.0]) rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -@test bugs_logp ≈ - LogDensityProblems.logdensity(t_p, [transform(bijector(dbin(0.1, 10)), 10.0)]) rtol = - 1E-6 diff --git a/test/logp_tests/single_distribution_models/ddirich.jl b/test/logp_tests/single_distribution_models/ddirich.jl deleted file mode 100644 index ae650c4fb..000000000 --- a/test/logp_tests/single_distribution_models/ddirich.jl +++ /dev/null @@ -1,46 +0,0 @@ -alpha = rand(10) -dist = ddirich(alpha) -test_θ_transformed = rand(9) -test_θ = DynamicPPL.invlink_and_reconstruct(dist, test_θ_transformed) - -rand(dist) -# Define the BUGS model -test_ddirich = @bugs begin - x[1:10] ~ ddirich(alpha[1:10]) -end - -# Compile the BUGS model -bugs_model = compile(test_ddirich, (alpha=alpha,), (x=test_θ,)) - -# Now, create a DynamicPPL model to represent the same distribution -@model function ddirich_test() - x ~ ddirich(alpha[1:10]) - return x -end - -dppl_model = ddirich_test() - -# Follow the rest of the structure of your script to test the log-density and other properties, -# similar to what you did for the LKJ distribution but adjusting for the Wishart distribution. - -p = DynamicPPL.LogDensityFunction(dppl_model) -t_p = DynamicPPL.LogDensityFunction( - dppl_model, - DynamicPPL.link!!(SimpleVarInfo(dppl_model), dppl_model), - DynamicPPL.DefaultContext(), -) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -dppl_logp = LogDensityProblems.logdensity(p, vcat(test_θ...)) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -bugs_logp_logp_ctx = evaluate!!( - JuliaBUGS.settrans(bugs_model, true), LogDensityContext(), test_θ_transformed -)[2] -@test bugs_logp ≈ bugs_logp_logp_ctx rtol = 1E-6 -@test bugs_logp ≈ LogDensityProblems.logdensity( - JuliaBUGS.settrans(bugs_model, true), test_θ_transformed -) rtol = 1E-6 -dppl_logp = LogDensityProblems.logdensity(t_p, test_θ_transformed) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/logp_tests/single_distribution_models/dwish.jl b/test/logp_tests/single_distribution_models/dwish.jl deleted file mode 100644 index ec0386d00..000000000 --- a/test/logp_tests/single_distribution_models/dwish.jl +++ /dev/null @@ -1,53 +0,0 @@ -# Define a scale matrix and degrees of freedom -scale_matrix = randn(10, 10) -scale_matrix = scale_matrix * transpose(scale_matrix) # Ensuring positive-definiteness -degrees_of_freedom = 12 - -dist = dwish(scale_matrix, degrees_of_freedom) -test_θ_transformed = rand(55) -test_θ = DynamicPPL.invlink_and_reconstruct(dist, test_θ_transformed) - -# Define the BUGS model -test_dwish = @bugs begin - x[1:10, 1:10] ~ dwish(scale_matrix[:, :], degrees_of_freedom) -end - -# Compile the BUGS model -bugs_model = compile( - test_dwish, - (degrees_of_freedom=degrees_of_freedom, scale_matrix=scale_matrix), - (x=test_θ,), -) - -# Now, create a DynamicPPL model to represent the same distribution -@model function dwish_test() - x ~ dwish(scale_matrix, degrees_of_freedom) - return x -end - -dppl_model = dwish_test() - -# Follow the rest of the structure of your script to test the log-density and other properties, -# similar to what you did for the LKJ distribution but adjusting for the Wishart distribution. - -p = DynamicPPL.LogDensityFunction(dppl_model) -t_p = DynamicPPL.LogDensityFunction( - dppl_model, - DynamicPPL.link!!(SimpleVarInfo(dppl_model), dppl_model), - DynamicPPL.DefaultContext(), -) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -dppl_logp = LogDensityProblems.logdensity(p, vcat(test_θ...)) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -bugs_logp_logp_ctx = evaluate!!( - JuliaBUGS.settrans(bugs_model, true), LogDensityContext(), test_θ_transformed -)[2] -@test bugs_logp ≈ bugs_logp_logp_ctx rtol = 1E-6 -@test bugs_logp ≈ LogDensityProblems.logdensity( - JuliaBUGS.settrans(bugs_model, true), test_θ_transformed -) rtol = 1E-6 -dppl_logp = LogDensityProblems.logdensity(t_p, test_θ_transformed) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/logp_tests/single_distribution_models/gamma.jl b/test/logp_tests/single_distribution_models/gamma.jl deleted file mode 100644 index 23e58512d..000000000 --- a/test/logp_tests/single_distribution_models/gamma.jl +++ /dev/null @@ -1,36 +0,0 @@ -model_def = @bugs begin - a ~ dgamma(0.001, 0.001) -end - -bugs_model = compile(model_def, NamedTuple(), (a=10,)) -vi = bugs_model.varinfo - -@model function dppl_gamma_model() - a ~ dgamma(0.001, 0.001) - return a -end - -dppl_model = dppl_gamma_model() - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -params_vi = JuliaBUGS.get_params_varinfo(bugs_model, vi) -# test if JuliaBUGS and DynamicPPL agree on parameters in the model -@test keys( - DynamicPPL.evaluate!!( - dppl_model, SimpleVarInfo(Dict{VarName,Any}()), DynamicPPL.SamplingContext() - )[2], -) == keys(params_vi) - -p = DynamicPPL.LogDensityFunction(dppl_model) -t_p = DynamicPPL.LogDensityFunction( - dppl_model, - DynamicPPL.link!!(SimpleVarInfo(dppl_model), dppl_model), - DynamicPPL.DefaultContext(), -) - -@test bugs_logp ≈ LogDensityProblems.logdensity(p, [10.0]) rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -@test bugs_logp ≈ - LogDensityProblems.logdensity(t_p, [transform(bijector(dgamma(0.001, 0.001)), 10.0)]) rtol = - 1E-6 diff --git a/test/logp_tests/single_distribution_models/lkj.jl b/test/logp_tests/single_distribution_models/lkj.jl deleted file mode 100644 index 03f2c6b1c..000000000 --- a/test/logp_tests/single_distribution_models/lkj.jl +++ /dev/null @@ -1,44 +0,0 @@ -dist = LKJ(10, 0.5) -test_θ_transformed = rand(45) -test_θ = DynamicPPL.invlink_and_reconstruct(dist, test_θ_transformed) - -test_lkj = @bugs begin - x[1:10, 1:10] ~ LKJ(10, 0.5) -end - -bugs_model = compile(test_lkj, NamedTuple(), (x=test_θ,)) -vi = bugs_model.varinfo - -# test param length given trans-dim bijectors -@test LogDensityProblems.dimension(JuliaBUGS.settrans(bugs_model, false)) == 100 -@test LogDensityProblems.dimension(JuliaBUGS.settrans(bugs_model, true)) == 45 - -@model function lkj_test() - x = Matrix{Float64}(undef, 10, 10) - x ~ LKJ(10, 0.5) - return x -end - -dppl_model = lkj_test() - -p = DynamicPPL.LogDensityFunction(dppl_model) -t_p = DynamicPPL.LogDensityFunction( - dppl_model, - DynamicPPL.link!!(SimpleVarInfo(dppl_model), dppl_model), - DynamicPPL.DefaultContext(), -) - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, false), DefaultContext())[2] -dppl_logp = LogDensityProblems.logdensity(p, vcat(test_θ...)) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 - -bugs_logp = JuliaBUGS.evaluate!!(JuliaBUGS.settrans(bugs_model, true), DefaultContext())[2] -bugs_logp_logp_ctx = evaluate!!( - JuliaBUGS.settrans(bugs_model, true), LogDensityContext(), test_θ_transformed -)[2] -@test bugs_logp ≈ bugs_logp_logp_ctx rtol = 1E-6 -@test bugs_logp ≈ LogDensityProblems.logdensity( - JuliaBUGS.settrans(bugs_model, true), test_θ_transformed -) rtol = 1E-6 -dppl_logp = LogDensityProblems.logdensity(t_p, test_θ_transformed) -@test bugs_logp ≈ dppl_logp rtol = 1E-6 diff --git a/test/runtests.jl b/test/runtests.jl index b62ce06b0..4bd3d8d2c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -53,7 +53,7 @@ if test_group == "compilation" || test_group == "all" include("bugs_primitives.jl") include("compile.jl") end - include("test_logp.jl") + include("log_density.jl") end if test_group == "gibbs" || test_group == "all" From ee7f64cc9a616dd7cf46453a0acd66fc4ca841e3 Mon Sep 17 00:00:00 2001 From: Xianda Sun Date: Mon, 7 Oct 2024 06:02:58 +0100 Subject: [PATCH 4/5] separate log_density test from compilation tests --- test/runtests.jl | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/test/runtests.jl b/test/runtests.jl index 4bd3d8d2c..7037d33e1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -23,7 +23,7 @@ using ReverseDiff AbstractMCMC.setprogress!(false) -const Tests = ("elementary", "compilation", "gibbs", "mcmchains", "all") +const Tests = ("elementary", "compilation", "log_density", "gibbs", "mcmchains", "all") const test_group = get(ENV, "TEST_GROUP", "all") if test_group ∉ Tests @@ -53,6 +53,9 @@ if test_group == "compilation" || test_group == "all" include("bugs_primitives.jl") include("compile.jl") end +end + +if test_group == "log_density" || test_group == "all" include("log_density.jl") end From 80fb74ca65f7831e5d9c1cdede1fe90a63cf33e5 Mon Sep 17 00:00:00 2001 From: Xianda Sun Date: Mon, 7 Oct 2024 06:20:08 +0100 Subject: [PATCH 5/5] remove extra imported names --- test/log_density.jl | 4 ---- 1 file changed, 4 deletions(-) diff --git a/test/log_density.jl b/test/log_density.jl index 3f1ed3489..9808208e0 100644 --- a/test/log_density.jl +++ b/test/log_density.jl @@ -1,7 +1,3 @@ -using JuliaBUGS: - BUGSGraph, DefaultContext, evaluate!!, get_params_varinfo, LogDensityContext -using DynamicPPL: DynamicPPL, getlogp, settrans!!, SimpleVarInfo, @model - # TODO: make this available in JuliaBUGS function _logjoint(model::JuliaBUGS.BUGSModel) return JuliaBUGS.evaluate!!(model, JuliaBUGS.DefaultContext())[2]