Factor Models
In this tutorial, we will cover a more advanced use-case of AdaptEllipticalSliceSampler.JL, where we use AGESS[1] to sample a subset of parameters, and use Gibbs sampling to conditionally update the rest of the samplers. In this tutorial, we consider the following factor model using the multiplicative gamma process shrinkage prior[2]:
\[\mathbf{y}_i = \Lambda \eta_i + \epsilon_i,\]
\[\epsilon_i \sim \mathcal{N}(0, \mathbf{D}),\]
\[\Lambda_{jh} \mid \phi_{jh}, \tau_h \sim \mathcal{N}(0, \phi_{jh}^{-1} \tau_h^{-1}),\]
\[\phi_{jh} \sim \text{Gamma}(\nu / 2, \nu / 2),\]
\[\tau_h = \prod_{k=1}^h \delta_k \;\;\; (h = 1, \dots, K)\]
\[\delta_1 \sim \text{Gamma}(a_1, 1),\]
\[\delta_{h} \sim \text{Gamma}(a_2, 1) \;\;\; (h \ge 2),\]
\[\mathbf{D}_{jj} \sim \text{Inv-Gamma}(a, b),\]
where $\mathbf{y}_i \in \mathbb{R}^P$, $\Lambda \in \mathbb{R}^{P \times K}$, $\eta_i \in \mathbb{R}^{K}$, $\mathbf{D}$ is a Diagonal matrix, and $a_2 > 1$ (to promote shrinkage). In these settings, it is often the case that $K$ is much smaller than $P$. Under this setup, a Gibbs sampler is available[2], making sampling relatively straightforward. Consider the case where we want to use AGESS to sample $\Lambda$, $\eta_i$, and $\mathbf{D}$, while using Gibbs updates for $\delta_h$ and $\phi_{jh}$. We will start by constructing Gibbs samplers for $\delta_h$ and $\phi_{jh}$.
Gibbs Updates
The conditional posterior distribution of $\delta_h$ and $\phi_{jh}$ are as follows[2]:
\[\phi_{jh}\mid \Theta_{-\phi_{jh}} \sim \text{Gamma}\left(0.5(\nu + 1), 0.5(\nu + \tau_h \Lambda_{jh}^2)\right)\]
\[\delta_1 \mid \Theta_{-\delta_1} \sim \text{Gamma}\left(a_1 + 0.5(P\times K), 1 + 0.5\left(\sum_{l=1}^K \tau_{l}^{(-1)} \sum_{j=1}^P \phi_{jl}\Lambda_{jl}^2\right)\right)\]
\[\delta_h \mid \Theta_{-\delta_h} \sim \text{Gamma}\left(a_2 + 0.5(P\times (K - h +1)), 1 + 0.5\left(\sum_{l=h}^K \tau_{l}^{(-h)} \sum_{j=1}^P \phi_{jl}\Lambda_{jl}^2\right)\right)\;\;\; (2 \le h \le K)\]
We can construct Julia functions to update these parameters as follows.
using Random, BenchmarkTools
using AdaptEllipticalSliceSampler
using Distributions
using Plots
using LinearAlgebra
Random.seed!(123)
function delta_sampler!(δ::AbstractVector{Y}, ϕ::AbstractMatrix{Y}, Λ::AbstractMatrix{Y},
a1_δ::Y, a2_δ::Y, ph_δ::AbstractVector{Y}) where {Y<:AbstractFloat}
K = length(δ)
P = size(Λ)[2]
## Sample δ1
α = a1_δ + 0.5 * K * P
@views ph_δ .= ϕ[1,:] .* Λ[1,:]
@views β = 1.0 + 0.5 * dot(ph_δ, Λ[1,:])
for i in 2:K
@views τ = prod(δ[2:i])
@views ph_δ .= ϕ[i,:] .* Λ[i,:]
@views β += 0.5 * τ * dot(ph_δ, Λ[i,:])
end
δ[1] = rand(Gamma(α, 1 / β))
## Sample rest of δ
for i in 2:K
α = a2_δ + 0.5 * (K - i + 1) * P
β = 1.0
@views τ = prod(δ[1:(i-1)])
for j in i:K
if i != j
τ *= δ[j]
end
@views ph_δ .= ϕ[j,:] .* Λ[j,:]
@views β += 0.5 * τ * dot(ph_δ, Λ[j,:])
end
δ[i] = rand(Gamma(α, 1 / β))
end
return nothing
end
function phi_sampler!(ϕ::AbstractMatrix{Y}, δ::AbstractVector{Y}, Λ::AbstractMatrix{Y},
ν::Y) where {Y<:AbstractFloat}
K = length(δ)
P = size(Λ)[2]
for i in 1:K
@views τ = prod(δ[1:i])
for j in 1:P
α = (ν + 1) * 0.5
β = 0.5 * (ν + τ * Λ[i,j]^2)
ϕ[i,j] = rand(Gamma(α, 1 / β))
end
end
return nothing
endphi_sampler! (generic function with 1 method)AGESS Updates
Now that we have our functions to update the $\phi$ and $\tau$ parameters via Gibbs updates, we can write a function that evaluates the conditional posterior of the rest of the parameters so that we can update the state of the Markov chain for these parameters via AGESS[1]. Here we will construct two functions. The first one evaluates the likelihood, and the second one reshapes and makes necessary transformations to our variables of interest.
function posterior(Λ::AbstractMatrix{Y}, η::AbstractMatrix{Y}, Y_obs::AbstractMatrix{Y},
D::AbstractMatrix{Y}, ϕ::AbstractMatrix{Y}, δ::AbstractVector{Y},
τ_ph::AbstractVector{Y}, a::Y, b::Y,
ph::AbstractVector{Y}, ph1::AbstractVector{Y})::Float64 where {Y<:AbstractFloat}
lpdf::Float64 = 0.0
## Likelihood
for i in 1:size(Y_obs)[1]
@views mul!(ph, Λ', η[i,:])
@views ph .-= Y_obs[i,:]
@views ph1 .= ph ./ D[diagind(D)]
@views lpdf += (-0.5* sum(log.(D[diagind(D)])) - 0.5 * dot(ph1, ph))
end
ig_d = InverseGamma(a,b)
for i in 1:size(D)[1]
lpdf += logpdf(ig_d, D[i,i])
end
τ_ph .= δ
for i in 2:length(δ)
τ_ph[i] *= τ_ph[i-1]
end
for i in 1:size(Λ)[1], j in 1:size(Λ)[2]
lpdf += 0.5 * log(ϕ[i,j] * τ_ph[i]) - 0.5 * (Λ[i,j]^2 * ϕ[i,j] * τ_ph[i])
end
for i in 1:size(η)[1]
@views lpdf += - 0.5 * dot(η[i,:], η[i,:])
end
return lpdf
end
function transform_posterior(x::AbstractVector{Y}, Y_obs::AbstractMatrix{Y},
ϕ::AbstractMatrix{Y}, δ::AbstractVector{Y},
τ_ph::AbstractVector{Y}, a::Y, b::Y,
N::T, P::T, K::T, D_ph::AbstractMatrix{Y},
ph::AbstractVector{Y},
ph1::AbstractVector{Y})::Float64 where {Y<:AbstractFloat, T<:Integer}
@views Λ_ph = reshape(x[1:K*P], (K, P))
@views η_ph = reshape(x[(K*P + 1):(N*K + K*P)], (N, K))
@views D_ph[diagind(D_ph)] .= exp.(x[(N*K + K*P + 1):(N*K + K*P + P)])
lpdf = posterior(Λ_ph, η_ph, Y_obs, D_ph, ϕ, δ, τ_ph, a, b, ph, ph1)
### Jacobian for transformation
@views lpdf += sum(x[(N*K + K*P + 1):(N*K + K*P + P)])
return lpdf
endtransform_posterior (generic function with 1 method)Construct Custom MCMC Sampling Scheme
Now that we have set up the functions necessary to use the AdaptEllipticalSliceSampler.jl package, and have constructed the functions necessary to perform Gibbs updates, we can use the following two functions to create a custom MCMC sampling method: AGESS_single_step_1d! and AGESS_single_step!. Here we will use a t-distribution with 6 degrees of freedom as our marginal distribution of the auxiliary variable in AGESS.
function custom_MCMC(Y_obs::AbstractMatrix{Y}, K::T, n_MCMC::T, a_1::Y, a_2::Y, a::Y, b::Y,
ϵ::Y, single_step_prop::Y, burnin::Y, ν::Y, β::Y) where {Y<:AbstractFloat, T<:Integer}
N = size(Y_obs)[1]
P = size(Y_obs)[2]
n_params = (N * K) + (K * P) + P
### Set up variables
x = zeros(n_MCMC, n_params)
ϕ = ones(n_MCMC, K, P)
δ = ones(n_MCMC, K)
## Allocate variables for AGESS
ph_AGESS = zeros(n_params)
z = similar(ph_AGESS)
w_const = max(2/3, ((cbrt(n_params) - 1) / cbrt(n_params)))
N_J = 2
n_j = 2
ph_cholesky_update = ones(n_params)
perm = randperm(n_params)
### Allocate variables for intermediate calculations
τ_ph = zeros(K)
D_ph = diagm(ones(P))
ph = zeros(P)
ph1 = similar(ph)
ph_δ = similar(ph)
### Setup adaptive parameters
μ_adapt = zeros(n_params)
μ_adapt_ph = zeros(n_params)
μ_0 = zeros(n_params)
Σ_adapt = diagm(ones(n_params))
Σ_chol_0 = cholesky(Σ_adapt)
Σ_chol_adapt = deepcopy(Σ_chol_0)
Σ_chol_adapt_ph = deepcopy(Σ_chol_0)
### variable for storing log posterior pdf
lpdf = zeros(n_MCMC)
@views lpdf[1] = transform_posterior(x[1,:], Y_obs, ϕ[1,:,:], δ[1,:], τ_ph, a,
b, N, P, K, D_ph, ph, ph1)
### Start MCMC
for i in 2:n_MCMC
### AGESS update
if i < (n_MCMC * burnin)
### 1-d AGESS updates for fast burn-in
@views lpdf[i] = AGESS_single_step_1d!(x, y -> transform_posterior(y, Y_obs,
ϕ[i,:,:], δ[i,:], τ_ph, a,
b, N, P, K, D_ph, ph, ph1), true, 6.0,
n_params, μ_adapt,
Σ_chol_adapt.L, lpdf[i-1], perm, i)
else
if rand() > (ϵ + single_step_prop)
### standard AGESS step
@views lpdf[i] = AGESS_single_step!(x, z, y -> transform_posterior(y, Y_obs,
ϕ[i,:,:], δ[i,:], τ_ph, a,
b, N, P, K, D_ph, ph, ph1), true, 6.0,
n_params, ph_AGESS, μ_adapt,
Σ_chol_adapt.L, lpdf[i-1], i)
elseif rand() < (single_step_prop / (ϵ + single_step_prop))
### AGESS updates with only 1-d updates
@views lpdf[i] = AGESS_single_step_1d!(x, y -> transform_posterior(y, Y_obs,
ϕ[i,:,:], δ[i,:], τ_ph, a,
b, N, P, K, D_ph, ph, ph1), true, 6.0,
n_params, μ_adapt,
Σ_chol_adapt.L, lpdf[i-1], perm, i)
else
### Non-adaptive update
@views lpdf[i] = AGESS_single_step!(x, z, y -> transform_posterior(y, Y_obs,
ϕ[i,:,:], δ[i,:], τ_ph, a,
b, N, P, K, D_ph, ph, ph1), true, 6.0,
n_params, ph_AGESS, μ_0,
Σ_chol_0.L, lpdf[i-1], i)
end
end
### Update variables
@views Λ_ph = reshape(x[i, 1:K*P], (K, P))
@views η_ph = reshape(x[(K*P + 1):(N*K + K*P)], (N, K))
@views D_ph[diagind(D_ph)] .= exp.(x[(N*K + K*P + 1):(N*K + K*P + P)])
### Gibbs Updates
@views phi_sampler!(ϕ[i,:,:], δ[i,:], Λ_ph, ν)
@views delta_sampler!(δ[i,:], ϕ[i,:,:], Λ_ph, a_1, a_2, ph_δ)
### Update Adaptive Scheme
w_i = i^(-w_const)
Σ_chol_adapt_ph.U .= sqrt((1 - w_i)) .* Σ_chol_adapt_ph.U
@views ph_cholesky_update .= sqrt(w_i) .* (x[i,:] .- μ_adapt_ph)
lowrankupdate!(Σ_chol_adapt_ph, ph_cholesky_update)
@views μ_adapt_ph .= (1 - w_i) * μ_adapt_ph + w_i * x[i,:]
## Adapt mean and covariance according to AIRMCMC
if i == N_J
Σ_chol_adapt.U .= Σ_chol_adapt_ph.U
@views μ_adapt .= μ_adapt_ph
n_j += 1
N_J += floor(n_j^β)
end
### Update next state
if i < n_MCMC
@views x[i+1,:] .= x[i,:]
@views ϕ[i+1,:,:] .= ϕ[i,:,:]
@views δ[i+1,:] .= δ[i,:]
end
### rerun evaluation of posterior density after gibbs updates
@views lpdf[i] = transform_posterior(x[i,:], Y_obs, ϕ[i,:,:], δ[i,:], τ_ph, a, b, N,
P, K, D_ph, ph, ph1)
### Print statement
if (i % 100) == 0
@views println("MCMC iter: ", i, " lpdf = ", mean(lpdf[i-99:i]))
end
end
Λ_out = zeros(n_MCMC, K, P)
η_out = zeros(n_MCMC, N, K)
D_out = zeros(n_MCMC, P, P)
for i in 1:n_MCMC
@views Λ_out[i,:,:] .= reshape(x[i, 1:K*P], (K, P))
@views η_out[i,:,:] .= reshape(x[i, (K*P + 1):(N*K + K*P)], (N, K))
for j in 1:P
D_out[i,j,j] = exp(x[i, (N*K + K*P + j)])
end
end
return Λ_out, η_out, D_out, ϕ, δ
endcustom_MCMC (generic function with 1 method)We have constructed our custom MCMC sampler. We can now generate some synthetic data and test the performance of our sampler.
### Generate Data
K = 4
N = 200
P = 50
Σ = diagm(ones(P))
ph = zeros(P)
for i in 1:(K - 3)
ph = randn(P)
Σ .+= i .* (ph * ph')
end
μ_0 = zeros(P)
Y_obs = rand(MvNormal(μ_0, Σ), N)'
Σ_truth = deepcopy(Σ)
### Run custom MCMC
a_1 = 2.0
a_2 = 2.0
ν = 2.0
a = 1.0
b = 1.0
n_MCMC = 10000
ϵ = 0.05
single_step_prop = 0.05
β = 0.5
burnin = 0.01
@time Λ_samp, η_samp, D_samp, ϕ_samp, δ_samp = custom_MCMC(Y_obs, K, n_MCMC, a_1, a_2, a,
b, ϵ, single_step_prop, burnin,
ν, β)([0.0 0.0 0.0 0.0; -0.11261588545665198 0.005249971151925579 0.16775604813056266 -0.2601908318136752; … ; -0.4783735333451412 0.14741277070885245 -0.13752651672205346 -0.7408121921811347; -0.4787869759209649 0.14775483731918365 -0.13837628298975976 -0.7413694684735557;;; 0.0 0.0 0.0 0.0; 0.6272518604464292 0.36343229606252503 0.06412071321084413 0.6000425306736885; … ; 0.661887397336742 0.1664980400663218 -0.19895723326924492 0.9332041479288784; 0.6626230992191375 0.16641821131918524 -0.19787325187540905 0.9329032366498556;;; 0.0 0.0 0.0 0.0; 0.0680130629945553 -0.11246160396480917 0.027213053442416955 0.22777614980501207; … ; 0.7986125595645264 0.14200420566155705 0.23203119105560258 0.7036986050720185; 0.7985216417599371 0.14322772904189673 0.23412380554844964 0.704969056428551;;; … ;;; 0.0 0.0 0.0 0.0; 0.02861140268179307 0.353416922916414 -0.06599431726184798 -0.4184013550528051; … ; 0.7823363026331739 0.44922533177954976 -0.38879821066788756 0.29242733480725563; 0.7839315580898109 0.44925757304307123 -0.38846200122125424 0.29226238326235826;;; 0.0 0.0 0.0 0.0; 0.8642754204290692 0.2744885939793009 0.10595661971511568 -0.2496774645556675; … ; 0.3906563216893849 0.8196303750413488 -0.05381090123950176 0.8577510986348574; 0.3931729610535239 0.8176001231396639 -0.053778016101519624 0.8573616861805332;;; 0.0 0.0 0.0 0.0; -0.17035259647226972 -0.6638863805273516 -0.1857441113331592 0.05076548849259288; … ; -0.7571764251174878 -0.1867554433250183 -0.244113463683035 -0.17581299917129034; -0.7543951551447285 -0.18817335778474778 -0.24380122838283672 -0.17492368350111023], [0.0 0.0 … 0.0 0.0; -0.3918981093771017 -0.13943981267369437 … 0.1379751099496585 0.5078759087300582; … ; -0.5019726974094655 -1.0288573247088004 … 0.67853122099909 0.7273370068068485; -0.49959219554990514 -1.0286293670643 … 0.6798159108354566 0.7275292080485456;;; 0.0 0.0 … 0.0 0.0; -0.466693301250591 -0.5980263234805476 … 1.2670680943774508 -0.0013491549191727691; … ; -0.006797541069831922 -0.5515830733992084 … 0.6755895546151077 0.15028687102243457; -0.007456751861762946 -0.5517624252017962 … 0.6788055781474722 0.14909319971316476;;; 0.0 0.0 … 0.0 0.0; -0.01421572823299255 0.04973088706271119 … -0.12319684678942582 -0.06214306411839359; … ; -0.011942301859877318 0.38523829539477333 … 0.36243503194815274 -0.13074559593477464; -0.01100963878069593 0.3850735180124295 … 0.35966306921823066 -0.1304763135327682;;; 0.0 0.0 … 0.0 0.0; -0.16390860705757668 -0.14852177707320424 … 0.008741964249883457 -0.2907388315004424; … ; -0.4304279455365867 -0.5743629592312818 … 0.8880590416844761 -0.5474500051676148; -0.4301927139338707 -0.5772901473467953 … 0.8888191350703712 -0.547378166009212], [1.0 0.0 … 0.0 0.0; 1.3138343896288092 0.0 … 0.0 0.0; … ; 1.1575286493071575 0.0 … 0.0 0.0; 1.15585894014701 0.0 … 0.0 0.0;;; 0.0 1.0 … 0.0 0.0; 0.0 1.0438001917752935 … 0.0 0.0; … ; 0.0 0.9956156831374872 … 0.0 0.0; 0.0 0.9953878484010158 … 0.0 0.0;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; … ;;; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0; … ; 0.0 0.0 … 0.0 0.0; 0.0 0.0 … 0.0 0.0;;; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 2.041801106997551 0.0; … ; 0.0 0.0 … 1.0622098627602783 0.0; 0.0 0.0 … 1.064897747479841 0.0;;; 0.0 0.0 … 0.0 1.0; 0.0 0.0 … 0.0 2.3308495889248704; … ; 0.0 0.0 … 0.0 1.0691863474064192; 0.0 0.0 … 0.0 1.0703406684455063], [1.0 1.0 1.0 1.0; 1.4538454644862804 1.4275743897916318 2.682847573938428 1.4883058539256195; … ; 0.5738550775992598 0.8815999100793224 0.15799864905536667 1.4253436781979612; 0.7694886587256188 2.3718776289186505 3.425807117732766 0.8572308838964827;;; 1.0 1.0 1.0 1.0; 1.2092593835693817 1.2936214327427986 0.9179504741479844 1.010729065181896; … ; 0.6582368563863007 2.575959239929949 0.4603431822028688 1.1600708074146937; 0.6774017625216484 0.11576996827076401 0.5566495847404886 0.13412892011848926;;; 1.0 1.0 1.0 1.0; 0.42463508672627126 3.0146916516943443 0.4904968513458115 1.1479821699368187; … ; 0.9516849979611975 0.8247757416567413 0.4426860415949487 0.10402745765323182; 1.1126850842598563 0.3660884236364071 0.6467145122680303 1.7025739794808508;;; … ;;; 1.0 1.0 1.0 1.0; 2.1831111037715294 1.1692777245747563 0.5319885506188426 1.1673298160364258; … ; 0.47730976978079226 0.9127279515027162 1.4657661299588982 1.8339145448154095; 2.762866231748155 1.372888897341869 0.45811020290972765 0.31045958991245826;;; 1.0 1.0 1.0 1.0; 1.0285802810213867 1.358044976682831 4.6766873822972785 0.35855623559499056; … ; 2.327986390013403 0.07018704874812852 0.25517142587012476 1.2294632644955084; 0.9481555919098233 0.34729988212393487 0.8374624392948276 0.42612432054629756;;; 1.0 1.0 1.0 1.0; 1.2990048864736647 1.206255639326305 0.13752134626669116 2.0222701427875616; … ; 0.39233742786068543 0.25433006968397087 1.145689960396191 0.33058201023562095; 0.6975721315635585 0.35345296396774806 0.5733256554528138 0.8409766302063405], [1.0 1.0 1.0 1.0; 4.882202470960154 1.5066751132266865 1.0132269540305536 0.8771241198653181; … ; 3.5738801534364346 3.2135532547317296 3.369343674305426 0.08961476451089931; 2.993787991654475 4.665606746495211 2.567963624895598 0.06538468818374196])We can look at estimates of $\Lambda'\Lambda + \mathbf{D}$ to see how well we recovered the covariance structure. First, we will look at the true covariance matrix.
heatmap(Σ_truth)
Now, we can look at the posterior element-wise mean of $\Lambda'\Lambda + \mathbf{D}$.
function posterior_Σ(Λ_samp::AbstractArray{Y,3}, D_samp::AbstractArray{Y,3}, P::T;
burnin = 0.5) where {Y<:AbstractFloat, T<:Integer}
n_MCMC = size(Λ_samp)[1]
burnin_num = floor(Int64, burnin * n_MCMC)
posterior_samps = zeros(n_MCMC - burnin_num, P, P)
for i in (burnin_num +1):n_MCMC
@views posterior_samps[i - burnin_num,:,:] .= Λ_samp[i,:,:]' * Λ_samp[i,:,:]
posterior_samps[i - burnin_num,:,:] .+= D_samp[i,:,:]
end
return posterior_samps
end
Σ_samp = posterior_Σ(Λ_samp, D_samp, P, burnin = 0.5)
Σ_mean = zeros(P,P)
for i in 1:P
for j in 1:P
Σ_mean[i,j] = mean(Σ_samp[:,i,j])
end
end
heatmap(Σ_mean)
We can also look at the trace plots of individual elements of $\Lambda'\Lambda + \mathbf{D}$, along with the true value (represented by the horizontal line).
plot(Σ_samp[:,1,1])
hline!([Σ_truth[1,1]])
Conclusion
As illustrated in this tutorial, the AdaptEllipticalSliceSampler.jl package can be used to construct custom MCMC schemes. When considering factor analysis in this setting, one may notice that we can simply use AGESS for sampling all parameters. Although this tutorial primarily serves as a guide for how to construct custom MCMC sampling schemes, we can compare the results obtained when using AGESS to sample all parameters.
function posterior2(Λ::AbstractMatrix{Y}, η::AbstractMatrix{Y}, Y_obs::AbstractMatrix{Y},
D::AbstractMatrix{Y}, ϕ::AbstractMatrix{Y}, δ::AbstractVector{Y},
τ_ph::AbstractVector{Y}, a_1::Y, a_2::Y, ν::Y, a::Y, b::Y,
ph::AbstractVector{Y}, ph1::AbstractVector{Y})::Float64 where {Y<:AbstractFloat}
lpdf::Float64 = 0.0
## Likelihood
for i in 1:size(Y_obs)[1]
@views mul!(ph, Λ', η[i,:])
@views ph .-= Y_obs[i,:]
@views ph1 .= ph ./ D[diagind(D)]
@views lpdf += (-0.5* sum(log.(D[diagind(D)])) - 0.5 * dot(ph1, ph))
end
##Priors
for i in eachindex(δ)
if i == 1
lpdf += logpdf(Gamma(a_1, 1), exp(δ[1])) + δ[1]
else
lpdf += logpdf(Gamma(a_2, 1), exp(δ[i])) + δ[i]
end
end
ig_d = InverseGamma(a,b)
for i in 1:size(D)[1]
lpdf += logpdf(ig_d, D[i,i]) + log(D[i,i])
end
τ_ph .= exp.(δ)
for i in 2:length(δ)
τ_ph[i] *= τ_ph[i-1]
end
gamma_d = Gamma(0.5 * ν, 2 / ν)
for i in 1:size(Λ)[1], j in 1:size(Λ)[2]
lpdf += 0.5 * log(exp(ϕ[i,j]) * τ_ph[i]) - 0.5 * (Λ[i,j]^2 * exp(ϕ[i,j]) * τ_ph[i])
lpdf += logpdf(gamma_d, exp(ϕ[i,j])) + ϕ[i,j]
end
for i in 1:size(η)[1]
@views lpdf += - 0.5 * dot(η[i,:], η[i,:])
end
return lpdf
end
function transform_posterior2(x::AbstractVector{Y}, Y_obs::AbstractMatrix{Y},
δ_ph::AbstractVector{Y}, τ_ph::AbstractVector{Y},
a_1::Y, a_2::Y, ν::Y, a::Y, b::Y,
N::T, P::T, K::T, D_ph::AbstractMatrix{Y},
ph::AbstractVector{Y},
ph1::AbstractVector{Y})::Float64 where {Y<:AbstractFloat, T<:Integer}
@views Λ_ph = reshape(x[1:K*P], (K, P))
@views η_ph = reshape(x[(K*P + 1):(N*K + K*P)], (N, K))
@views D_ph[diagind(D_ph)] .= exp.(x[(N*K + K*P + 1):(N*K + K*P + P)])
@views ϕ_ph = reshape(x[(N*K + K*P + P + 1):(N*K + 2*K*P + P)], (K, P))
@views δ_ph .= x[(N*K + 2*K*P + P + 1):(N*K + 2*K*P + P + K)]
lpdf = posterior2(Λ_ph, η_ph, Y_obs, D_ph, ϕ_ph, δ_ph, τ_ph, a_1, a_2, ν, a, b, ph, ph1)
return lpdf
endtransform_posterior2 (generic function with 1 method)Now that we have specified a function evaluating the posterior log pdf, we can call the function AGESS.
# Set variables
τ_ph = zeros(K)
δ_ph = zeros(K)
D_ph = zeros(P,P)
ph = zeros(P)
ph1 = similar(ph)
n_params = N*K + 2*K*P + P + K
results = AGESS(x -> transform_posterior2(x, Y_obs, δ_ph, τ_ph, a_1, a_2, ν, a, b, N, P, K,
D_ph, ph, ph1), n_MCMC, n_params)AdaptEllipticalSliceSampler.MCMC_output([0.0 0.0 … 0.0 0.0; -0.43065264858604724 0.12389865198690232 … 0.4034177716876439 -0.15914288288687442; … ; -0.3047840834660331 -0.4373033921721066 … 0.38168489403328737 -0.8767973544541903; -0.304889413474845 -0.43747364882834505 … 0.38165380466229637 -0.8767515985795034], [-9792.603397906596, -7726.336233723766, -6849.226703702453, -6274.567580506465, -5934.757641962126, -5832.48164667897, -5747.620801721173, -5698.669186667385, -5609.2753344062485, -5561.010747446548 … -4998.009021008172, -4998.317680309515, -4999.013509905073, -4998.861432427295, -4998.60154080151, -4999.199104639747, -4998.92005593151, -4998.681408026886, -4998.246179067437, -4998.151222782084], AGESS_MCMC_params(Main.var"#7#8"(), [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0], [1.0 0.0 … 0.0 0.0; 0.0 1.0 … 0.0 0.0; … ; 0.0 0.0 … 1.0 0.0; 0.0 0.0 … 0.0 1.0], true, 6.0, 0.5, 0.05, 0.5, 0.05, 10000, 1254), [0.0025366550824784994 -9.838165253706409e-5 … -6.482639483650268e-6 -5.1028471232458524e-5; -9.838165253706409e-5 0.0016796865602280095 … 9.243397589921745e-6 -0.00012378207487112318; … ; -6.482639483650268e-6 9.243397589921745e-6 … 0.0004867475201965747 -4.835504777845273e-5; -5.1028471232458524e-5 -0.00012378207487112318 … -4.8355047778452724e-5 0.0005016575776250623], [-0.3696602193126292, -0.3156574718927191, -0.39862928761166244, -0.3760652143364195, 0.6524130626661507, 0.44222381879412914, 0.6020073016986791, 0.3970015136457124, 0.4750600096110217, 0.5489646167425547 … 0.24540710373202085, -0.9498564245473958, -0.4916184851081164, 0.03587274254842831, 0.01820905549740163, 0.05895132597165255, 1.197738908984954, 0.9407767100753788, 0.3662812400292673, -0.9214111676815309])We can extract the parameters using the following code.
Λ_samp2 = zeros(n_MCMC, K, P)
η_samp2 = zeros(n_MCMC, N, K)
D_samp2 = zeros(n_MCMC, P, P)
δ_samp2 = zeros(n_MCMC, K)
ϕ_samp2 = zeros(n_MCMC, K, P)
for i in 1:n_MCMC
@views Λ_samp2[i,:,:] .= reshape(results.samps[i, 1:K*P], (K, P))
@views η_samp2[i,:,:] .= reshape(results.samps[i, (K*P + 1):(N*K + K*P)], (N, K))
for j in 1:P
D_samp2[i,j,j] = exp(results.samps[i, (N*K + K*P + j)])
end
@views ϕ_samp2[i,:,:] .= reshape(exp.(results.samps[i, (N*K + K*P + P + 1):(N*K + 2*K*P + P)]), (K, P))
@views δ_samp2[i,:] .= exp.(results.samps[i, (N*K + 2*K*P + P + 1):(N*K + 2*K*P + P + K)])
endSimilarly to before, we can look at the posterior element-wise mean of $\Lambda'\Lambda + \mathbf{D}$.
Σ_samp2 = posterior_Σ(Λ_samp2, D_samp2, P, burnin = 0.5)
Σ_mean2 = zeros(P,P)
for i in 1:P
for j in 1:P
Σ_mean2[i,j] = mean(Σ_samp2[:,i,j])
end
end
heatmap(Σ_mean2)
Lastly, we can look at the trace plots of individual elements of $\Lambda'\Lambda + \mathbf{D}$, along with the true value (represented by the horizontal line).
plot(Σ_samp2[:,1,1])
hline!([Σ_truth[1,1]])