made sure that pgas can sample t=0 properly. Updated example scripts

Author: mattiasvillani

.github/workflows/CI.yml

@@ -24,7 +24,6 @@ jobs:
       matrix:
         version:
           - '1.11'
-          - 'pre'
         os:
           - ubuntu-latest
         arch:

docs/src/examples/LGSS.md

@@ -27,7 +27,6 @@ gr(legend = :topleft, grid = false, color = colors[2], lw = 2, legendfontsize=8,
     xtickfontsize=8, ytickfontsize=8, xguidefontsize=8, yguidefontsize=8,
     titlefontsize = 10, markerstrokecolor = :auto)
 
-myquantile(A, p; dims, kwargs...) = mapslices(x -> quantile(x, p; kwargs...), A; dims)
 Random.seed!(123);
 
 # ### Set up the state-space model
@@ -38,7 +37,7 @@ mutable struct LGSSParams
     σₑ::Float64
 end
 
-prior(θ) = Normal(0, θ.σᵥ / √((1 - θ.a^2)));
+prior(θ) = Normal(0, 10*θ.σᵥ / √((1 - θ.a^2)));
 transition(θ, state, t) = Normal(θ.a * state, θ.σᵥ);  
 observation(θ, state, t) = Normal(state, θ.σₑ);   
 
@@ -65,25 +64,26 @@ plot!(y; seriestype=:scatter, label="observed, y", xlabel="t", markersize = 2,
     color = colors[1], markerstrokecolor = :auto)
 
 # ### PGAS sampling
-Nₚ = 20      # Number of particles for PGAS
-Nₛ = 1000;   # Number of samples from posterior
-PGASdraws = PGASsampler(y, θ, Nₛ, Nₚ, prior, transition, observation)
+Nₚ = 20             # Number of particles for PGAS
+Nₛ = 10000;         # Number of samples from posterior
+sample_t0 = true    # Sample state at t=0 ?
+PGASdraws = PGASsampler(y, θ, Nₛ, Nₚ, prior, transition, observation; sample_t0 = sample_t0)
 PGASmean = mean(PGASdraws, dims = 3)[:,:,1]
-PGASquantiles = myquantile(PGASdraws, [0.025, 0.975], dims = 3);
+PGASquantiles = quantile_multidim(PGASdraws, [0.025, 0.975], dims = 3);
 
 # ### FFBS Sampling
 ## Set up the LGSS for FFBS and sample
 Σₑ = [θ.σₑ^2]
 Σₙ = [θ.σᵥ^2]
 μ₀ = [0;;]
-Σ₀ = [θ.σᵥ^2/(1-θ.a^2);;]
+Σ₀ = 10^2*[θ.σᵥ^2/(1-θ.a^2);;]
 A = θ.a
 C = 1
 B = 0
 U = zeros(T,1);
-FFBSdraws = FFBS(U, y, A, B, C, Σₑ, Σₙ, μ₀, Σ₀, Nₛ);
-FFBSmean = mean(FFBSdraws, dims = 3)[2:end,:,1] # Exclude initial state at t=0
-FFBSquantiles = myquantile(FFBSdraws, [0.025, 0.975], dims = 3)[2:end,:,:];
+FFBSdraws = FFBS(U, y, A, B, C, Σₑ, Σₙ, μ₀, Σ₀, Nₛ; sample_t0 = sample_t0);
+FFBSmean = mean(FFBSdraws, dims = 3) 
+FFBSquantiles = quantile_multidim(FFBSdraws, [0.025, 0.975], dims = 3);
 
 # ### Plot the posterior mean and 95% C.I. intervals from both algorithms
 plottrue = true
@@ -93,7 +93,7 @@ for j in 1:p
 
     #True state evolution
     if plottrue
-        plt_tmp = plot(x, c = colors[3], lw = 1, label = "true state")
+        plt_tmp = plot([NaN*ones(sample_t0);x], c = colors[3], lw = 1, label = "true state")
     else
         plt_tmp = plot()
     end

docs/src/examples/Poisson.md

@@ -24,7 +24,6 @@ gr(legend = :topleft, grid = false, color = colors[2], lw = 2, legendfontsize=8,
     xtickfontsize=8, ytickfontsize=8, xguidefontsize=8, yguidefontsize=8,
     titlefontsize = 10, markerstrokecolor = :auto)
 
-myquantile(A, p; dims, kwargs...) = mapslices(x -> quantile(x, p; kwargs...), A; dims)
 Random.seed!(123);
 
 # ### Set up Poisson model
@@ -79,7 +78,7 @@ function plotEvolDistributions!(plt, postDraws, quantiles, trueEvol = nothing,
         label = nothing, shaded = false; plotSettings...)
 
     T, nState, nSim = size(postDraws)
-    postquantiles = myquantile(postDraws, quantiles, dims = 3);
+    postquantiles = quantile_multidim(postDraws, quantiles, dims = 3);
 
     if !isnothing(trueEvol)
         plot!(plt, 1:T, trueEvol, lw = 1, c = :black, label = "true"; plotSettings...)
@@ -122,17 +121,20 @@ plt
 # ### PGAS sampling
 nSim = 1000;             #  Number of samples from posterior
 nParticles = 100         # Number of particles for PGAS
-PGASdraws = PGASsampler(y, θ, nSim, nParticles, prior, transition, observation);
+sample_t0 = true         # Sample state at t=0 ?
+PGASdraws = PGASsampler(y, θ, nSim, nParticles, prior, transition, observation; 
+    sample_t0 = sample_t0);
 
 # Plot the true evolution and the posterior distributions from PGAS
 quantiles = [0.025, 0.5, 0.975]
 pltLogIntensity = plot(; title = "Log Intensity "*L" \log\lambda_t = x_t")
-plotEvolDistributions!(pltLogIntensity, PGASdraws, quantiles, x, 
+plotEvolDistributions!(pltLogIntensity, PGASdraws, quantiles, [NaN*ones(sample_t0);x], 
     "PGAS(N=$nParticles)", true; color = colors[3], lw = 1, legend = :topright)
 
 pltIntensity = plot(; title = "Intensity "*L" \lambda_t = \exp(x_t)")
-plotEvolDistributions!(pltIntensity, exp.(PGASdraws), quantiles, exp.(x), 
-    "PGAS(N=$nParticles)", true; color = colors[3], lw = 1, legend = :topright)
+plotEvolDistributions!(pltIntensity, exp.(PGASdraws), quantiles, 
+    exp.([NaN*ones(sample_t0);x]), "PGAS(N=$nParticles)", true; color = colors[3], 
+    lw = 1, legend = :topright)
 
 plot(pltLogIntensity, pltIntensity, layout = (1,2), size = (800, 300), bottommargin = 5mm)
 
@@ -147,7 +149,8 @@ B = θ.μ*(1-θ.a) # The transition model is x_t = μ + a (x_{t-1} - μ) + η_t,
 U = ones(T,1);
 
 # Simulate from the Laplace approximation of the posterior distribution 
-LaplaceDraws = FFBS_laplace(U, y, A, B, Σₙ, μ₀, Σ₀, observation, θ, nSim);
+LaplaceDraws = FFBS_laplace(U, y, A, B, Σₙ, μ₀, Σ₀, observation, θ, nSim; 
+    sample_t0 = sample_t0);
 
 # Plot the PGAS and Laplace posterior distributions, this time without true evolution
 pltLogIntensity = plot(; title = "Log Intensity "*L" \log\lambda_t = x_t")
@@ -179,10 +182,12 @@ x, y, plt = simTvPoisson(observation, transition, prior, θ, T);
 plt
 
 # Simulate from the PGAS posterior distribution
-PGASdraws = PGASsampler(y, θ, nSim, nParticles, prior, transition, observation);
+PGASdraws = PGASsampler(y, θ, nSim, nParticles, prior, transition, observation; 
+    sample_t0 = sample_t0);
 
 # Simulate from the Laplace approximation of the posterior distribution 
-LaplaceDraws = FFBS_laplace(U, y, A, B, Σₙ, μ₀, Σ₀, observation, θ, nSim);
+LaplaceDraws = FFBS_laplace(U, y, A, B, Σₙ, μ₀, Σ₀, observation, θ, nSim; 
+    sample_t0 = sample_t0);
 
 # Plot the PGAS and Laplace posterior distributions, without true evolution
 pltLogIntensity = plot(; title = "Log Intensity "*L" \log\lambda_t = x_t")

docs/src/examples/SAR.md

@@ -2,7 +2,7 @@
 
 julia_version = "1.11.6"
 manifest_format = "2.0"
-project_hash = "69971822c41d250d3341c5f32e3985a0a89a97bb"
+project_hash = "1442dc65b9337561130506a3e2dc47606ca8686b"
 
 [[deps.ADTypes]]
 git-tree-sha1 = "7927b9af540ee964cc5d1b73293f1eb0b761a3a1"
@@ -297,11 +297,6 @@ version = "0.7.4"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
     Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
-[[deps.DisplayAs]]
-git-tree-sha1 = "43c017d5dd3a48d56486055973f443f8a39bb6d9"
-uuid = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6"
-version = "0.1.6"
-
 [[deps.Distributed]]
 deps = ["Random", "Serialization", "Sockets"]
 uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"