From e5c578c422dee7c9fcd3b100abe326e8abb64f57 Mon Sep 17 00:00:00 2001 From: mschauer Date: Thu, 25 Aug 2022 17:14:37 +0200 Subject: [PATCH 1/3] Bouncy particle sampler with ZigZagBoomerang.jl --- Julia/Project.toml | 9 ++++++ Julia/bouncy.jl | 75 ++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 84 insertions(+) create mode 100644 Julia/Project.toml create mode 100644 Julia/bouncy.jl diff --git a/Julia/Project.toml b/Julia/Project.toml new file mode 100644 index 0000000..f7c80c9 --- /dev/null +++ b/Julia/Project.toml @@ -0,0 +1,9 @@ +[deps] +DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +ParquetFiles = "46a55296-af5a-53b0-aaa0-97023b66127f" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" +ZigZagBoomerang = "36347407-b186-4a6a-8c98-4f4567861712" diff --git a/Julia/bouncy.jl b/Julia/bouncy.jl new file mode 100644 index 0000000..86c4c08 --- /dev/null +++ b/Julia/bouncy.jl @@ -0,0 +1,75 @@ +# Sample a logistic regression poster (Pima data) using the bouncy particle sampler +# with https://github.com/mschauer/ZigZagBoomerang.jl +using ZigZagBoomerang +using ParquetFiles, DataFrames, Random, Distributions, StatsBase, ForwardDiff, Plots +using LinearAlgebra + +# Load and process the data +df = DataFrame(load(joinpath(@__DIR__, "..", "pima.parquet"))) +y = map((yi) -> yi == "Yes" ? 1.0 : 0.0, df.type) +X = hcat(ones(200, 1), Matrix(df[:,1:7])) + +# Define target density/potential +function lprior(β) + logpdf(Normal(0,10), β[1]) + sum(bi -> logpdf(Normal(0,1), bi), @view β[2:8]) +end +ll(β, X, y) = sum(-log.(1.0 .+ exp.(-(2*y .- 1.0).*(X * β)))) +ϕ(β, X, y) = -(lprior(β) + ll(β, X, y)) # potential + +# Actually, we need a gradient... +function ∇ϕ!(out, t, β, _, _, x, y) # gradient of ϕ + @. out = -β + out[1] = out[1]/100 # scale + out .+= transpose(x) * (y .- (1.0 ./ (1.0 .+ exp.(-x * β)))) + out +end + +# ... and directional derivatives we get with ForwardDiff +function dϕ(t, x, v, _, X, y) # two directional derivatives of ϕ + u = ForwardDiff.derivative(t -> ϕ(x + t*v, X, y), ForwardDiff.Dual{:dϕ}(0.0, 1.0)) + u.value, u.partials[] +end + +# Diagonal covariance estimate from warm-up (en lieu of mass matrix) +M = ZigZagBoomerang.PDMats.PDiagMat([3.5, 0.0038, 4.4e-5, 0.00036, 0.00052, 0.002, 0.33, 0.00048]) + +# Going to use BouncyParticle sampler + +Z = BouncyParticle(missing, missing, # ignored + 2.0, # momentum refreshment rate + 0.9, # momentum correlation / only gradually change momentum in refreshment/momentum update + M, # metric + missing +) + +# Initialize sampler +d = 8 # number of parameters +x0 = zeros(d) # starting point sampler +θ0 = randn(d) # starting direction sampler +T = 4000. # end time (similar to number of samples in MCMC) +bound = ZigZagBoomerang.LocalBound(20.0) # local bound on the rate of direction changes + + +# Call the sampler +trace, final, (acc, num), cs = @time pdmp( + dϕ, # return first two directional derivatives of negative target log-density in direction v + ∇ϕ!, # return gradient of negative target log-density + 0.0, x0, θ0, T, # initial state and duration + bound, # inital guess for bound + Z, # sampler + X, y; # data + progress=true, # show progress bar +) +# Get time and location of sampling points seperately as two vectors +ts, xs = ZigZagBoomerang.sep(trace) + +# Plots.jl wants matrix and not vector of vectors +out = [x[j] for x in xs, j in 1:d] + +# Plot results +plot(ts, out, layout=(4, 2)) +savefig("tracebouncy.pdf") +histogram(out, layout=(4, 2)) +savefig("histbouncy.pdf") +plot(1:400, autocor(out, 1:400), layout = (4, 2)) +savefig("acfbouncy.pdf") \ No newline at end of file From 53e54332cdf8e72c43ae41e7e114c452cd598e32 Mon Sep 17 00:00:00 2001 From: mschauer Date: Fri, 26 Aug 2022 12:24:47 +0200 Subject: [PATCH 2/3] Cosmetics --- Julia/bouncy.jl | 40 ++++++++++++++++++++-------------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/Julia/bouncy.jl b/Julia/bouncy.jl index 86c4c08..7a056b3 100644 --- a/Julia/bouncy.jl +++ b/Julia/bouncy.jl @@ -9,24 +9,23 @@ df = DataFrame(load(joinpath(@__DIR__, "..", "pima.parquet"))) y = map((yi) -> yi == "Yes" ? 1.0 : 0.0, df.type) X = hcat(ones(200, 1), Matrix(df[:,1:7])) +# Prior +prior = MvNormal(Diagonal([10.0; ones(7)].^2)) + # Define target density/potential -function lprior(β) - logpdf(Normal(0,10), β[1]) + sum(bi -> logpdf(Normal(0,1), bi), @view β[2:8]) -end ll(β, X, y) = sum(-log.(1.0 .+ exp.(-(2*y .- 1.0).*(X * β)))) -ϕ(β, X, y) = -(lprior(β) + ll(β, X, y)) # potential +ϕ(β, prior, X, y) = -(logpdf(prior, β) + ll(β, X, y)) # potential # Actually, we need a gradient... -function ∇ϕ!(out, t, β, _, _, x, y) # gradient of ϕ - @. out = -β - out[1] = out[1]/100 # scale - out .+= transpose(x) * (y .- (1.0 ./ (1.0 .+ exp.(-x * β)))) +function ∇ϕ!(out, t, β, _, _, prior, X, y) # gradient of ϕ + @. out = -β/prior.Σ.diag + out .+= transpose(X) * (y .- (1 ./ (1 .+ exp.(-X * β)))) out end # ... and directional derivatives we get with ForwardDiff -function dϕ(t, x, v, _, X, y) # two directional derivatives of ϕ - u = ForwardDiff.derivative(t -> ϕ(x + t*v, X, y), ForwardDiff.Dual{:dϕ}(0.0, 1.0)) +function dϕ(t, x, v, _, prior, X, y) # two directional derivatives of ϕ + u = ForwardDiff.derivative(t -> ϕ(x + t*v, prior, X, y), ForwardDiff.Dual{:dϕ}(0.0, 1.0)) u.value, u.partials[] end @@ -34,7 +33,6 @@ end M = ZigZagBoomerang.PDMats.PDiagMat([3.5, 0.0038, 4.4e-5, 0.00036, 0.00052, 0.002, 0.33, 0.00048]) # Going to use BouncyParticle sampler - Z = BouncyParticle(missing, missing, # ignored 2.0, # momentum refreshment rate 0.9, # momentum correlation / only gradually change momentum in refreshment/momentum update @@ -49,7 +47,6 @@ x0 = zeros(d) # starting point sampler T = 4000. # end time (similar to number of samples in MCMC) bound = ZigZagBoomerang.LocalBound(20.0) # local bound on the rate of direction changes - # Call the sampler trace, final, (acc, num), cs = @time pdmp( dϕ, # return first two directional derivatives of negative target log-density in direction v @@ -57,19 +54,22 @@ trace, final, (acc, num), cs = @time pdmp( 0.0, x0, θ0, T, # initial state and duration bound, # inital guess for bound Z, # sampler - X, y; # data + prior, X, y; # data progress=true, # show progress bar ) -# Get time and location of sampling points seperately as two vectors -ts, xs = ZigZagBoomerang.sep(trace) +# Time and location seperately as two vectors (instead of of a vector of pairs) +ts, xs = ZigZagBoomerang.splitpairs(trace) -# Plots.jl wants matrix and not vector of vectors +# Plots.jl wants a matrix and not vector of vectors out = [x[j] for x in xs, j in 1:d] # Plot results -plot(ts, out, layout=(4, 2)) +labels = ["β$i" for _ in 1:1, i in 1:d] +plot(ts, out, layout=(4, 2), label=labels) savefig("tracebouncy.pdf") -histogram(out, layout=(4, 2)) +histogram(out, layout=(4, 2), label=labels) savefig("histbouncy.pdf") -plot(1:400, autocor(out, 1:400), layout = (4, 2)) -savefig("acfbouncy.pdf") \ No newline at end of file +plot(1:200, autocor(out, 1:200), label=labels, layout = (4, 2)) +savefig("acfbouncy.pdf") +plot(getindex.(xs,1), getindex.(xs,2), xlabel= labels[1], ylabel=labels[2], label="trace") +savefig("pairbouncy.pdf") \ No newline at end of file From 7d7404f2057bd1d673f39cd2d841be98e821ca09 Mon Sep 17 00:00:00 2001 From: mschauer Date: Fri, 26 Aug 2022 12:36:54 +0200 Subject: [PATCH 3/3] tweaks --- Julia/bouncy.jl | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/Julia/bouncy.jl b/Julia/bouncy.jl index 7a056b3..05e2c94 100644 --- a/Julia/bouncy.jl +++ b/Julia/bouncy.jl @@ -7,7 +7,8 @@ using LinearAlgebra # Load and process the data df = DataFrame(load(joinpath(@__DIR__, "..", "pima.parquet"))) y = map((yi) -> yi == "Yes" ? 1.0 : 0.0, df.type) -X = hcat(ones(200, 1), Matrix(df[:,1:7])) +d = 8 # number of parameters +X = hcat(ones(length(y), 1), Matrix(df[:,1:d-1])) # Prior prior = MvNormal(Diagonal([10.0; ones(7)].^2)) @@ -17,15 +18,15 @@ ll(β, X, y) = sum(-log.(1.0 .+ exp.(-(2*y .- 1.0).*(X * β)))) ϕ(β, prior, X, y) = -(logpdf(prior, β) + ll(β, X, y)) # potential # Actually, we need a gradient... -function ∇ϕ!(out, t, β, _, _, prior, X, y) # gradient of ϕ +function ∇ϕ!(out, t, β, v, _, prior, X, y) # gradient of ϕ @. out = -β/prior.Σ.diag out .+= transpose(X) * (y .- (1 ./ (1 .+ exp.(-X * β)))) out end # ... and directional derivatives we get with ForwardDiff -function dϕ(t, x, v, _, prior, X, y) # two directional derivatives of ϕ - u = ForwardDiff.derivative(t -> ϕ(x + t*v, prior, X, y), ForwardDiff.Dual{:dϕ}(0.0, 1.0)) +function dϕ(t, β, v, _, prior, X, y) # two directional derivatives of ϕ + u = ForwardDiff.derivative(t -> ϕ(β + t*v, prior, X, y), ForwardDiff.Dual{:dϕ}(0.0, 1.0)) u.value, u.partials[] end @@ -34,14 +35,13 @@ M = ZigZagBoomerang.PDMats.PDiagMat([3.5, 0.0038, 4.4e-5, 0.00036, 0.00052, 0.00 # Going to use BouncyParticle sampler Z = BouncyParticle(missing, missing, # ignored - 2.0, # momentum refreshment rate - 0.9, # momentum correlation / only gradually change momentum in refreshment/momentum update + 0.5, # momentum refreshment rate + 0.7, # momentum correlation / only gradually change momentum in refreshment/momentum update M, # metric missing ) # Initialize sampler -d = 8 # number of parameters x0 = zeros(d) # starting point sampler θ0 = randn(d) # starting direction sampler T = 4000. # end time (similar to number of samples in MCMC) @@ -69,7 +69,7 @@ plot(ts, out, layout=(4, 2), label=labels) savefig("tracebouncy.pdf") histogram(out, layout=(4, 2), label=labels) savefig("histbouncy.pdf") -plot(1:200, autocor(out, 1:200), label=labels, layout = (4, 2)) +scatter(1:100, autocor(out, 1:100), label=labels, layout = (4, 2)) savefig("acfbouncy.pdf") plot(getindex.(xs,1), getindex.(xs,2), xlabel= labels[1], ylabel=labels[2], label="trace") savefig("pairbouncy.pdf") \ No newline at end of file