0
$\begingroup$

In the Pseudo-marginal Metropolis-Hastings algorithm exact sampling of a posterior distribution is performed when using an unbiased estimate of the marginal likelihood. However, I am having problems with the unbiasedness of the marginal likelihood.

Consider the model $$ \mathbf{y}_n = \begin{bmatrix} 1 \\ 1\end{bmatrix}s_n + \begin{bmatrix} 0 \\ 1\end{bmatrix} a + \mathbf{e}_n $$ where $\mathbf{e}_n \sim \mathcal{N}(0, \Sigma_e)$. $a$ is the parameter of interest and $s_n$ is the nuisance parameter. For $N$ measurements we have $$ \mathbf{y} = \mathbf{g} a + \mathbf{H} \mathbf{s}+ \mathbf{e}, $$ where $$\mathbf{g} = 1_N \otimes \begin{bmatrix} 0 \\ 1\end{bmatrix} \quad \mathbf{H} = I_N \otimes \begin{bmatrix} 1 \\ 1\end{bmatrix} $$ and $\mathbf{y}$ are the stacked measurements and $\mathbf{e}$ and $\mathbf{s}$ similarly defined. Note that $\mathbf{e} \sim \mathcal{N}(\mathbf{0}, \mathbf{Q})$ where $\mathbf{Q} = I_N \otimes \Sigma_e$.

Assume the prior of $\mathbf{s}$ is $\mathcal{N}( \mathbf{s} | \mathbf{\mu}_s, \Sigma_s)$. The marginal likelihood given the parameter $a$ $$ p(\mathbf{y}|a) = \int p(\mathbf{y}, \mathbf{s}|a) d \mathbf{s} = \mathcal{N}(\mathbf{y}| \mathbf{g}a + \mathbf{H}\mathbf{\mu}_s, \mathbf{Q} + \mathbf{H} \Sigma_s \mathbf{H}^\top) $$ and the pseudo marginal likelihood is (using $p(\mathbf{y}, \mathbf{s}|a) = p(\mathbf{y}| \mathbf{s},a) p(\mathbf{s} | a)$) $$ \widehat{p(\mathbf{y}|a)} = \frac{1}{M} \sum_{m=1}^{M} p(\mathbf{y}| \mathbf{s}^{(m)},a) \\ \mathbf{s}^{(m)} \sim p(\mathbf{s} ) = \mathcal{N}( \mathbf{s} | \mathbf{\mu}_s, \Sigma_s) $$ So I expected the pseudo marginal likelihood $\widehat{p(\mathbf{y}|a)}$ being an unbiased estimate of $p(\mathbf{y}|a)$ for all $M$, and the variance of $\widehat{p(\mathbf{y}|a)}$ will depend on $M$.

However, when I try to simulate I get a bias for low $M$, Mean and variance The distribution of the estimates is shown below distribution

The code is given below:

using Distributions
using LinearAlgebra
using PyPlot
using Printf

pygui(true)

eye(N) =  Matrix{Float64}(I,N,N)
⊗(A,B) = kron(A,B)

N = 10
g = ones(N) ⊗ [0.0; 1.0]
H = eye(N) ⊗ ones(2)

a = 1.0
S_prior = MvNormal(zeros(N), 2)
s = rand(S_prior)
Y = MvNormal(g*a + H*s,1)
y = rand(Y)
a0 = 1.01

p_y_a(a) = MvNormal(g*a + H*mean(S_prior), cov(Y) + H*cov(S_prior)*H')
p_y_a_samples(a, M) = begin
    map(1:M) do m
        s_m = rand(S_prior)
        p_y_s_a = MvNormal(g*a + H*s_m, cov(Y))
        pdf(p_y_s_a, y)
    end
end

M_tot = [1, 2, 4, 10, 100, 1000, 10_000]
M_samples = map(n ->  Int(1e5/n),N_tot)
samples = map(M_samples) do m
    Array{Float64}(undef, m)
end
let
    @progress for (k,M) in enumerate(M_tot)
        for m = 1:N_samples[k]
            samples[k][m] = log(mean(p_y_a_samples(a0, M)))
        end
    end
end
fig, ax = subplots(1,1)
vars = map(var, samples)
for k = 1:length(M_tot)
    ax.hist(samples[k], bins = "auto", histtype = "step",
            label = @sprintf(
                "Np = %d, # realizations = %d,var = %3.f", M_tot[k], N_samples[k], vars[k]
            )
    )
end
ax.grid(true)
ax.legend(loc = "upper left")
ax.set_xlabel("log(p_hat(y|a))")
ax.set_ylabel("Count")

fig, ax = subplots(2,1)
means = map(mean, samples)
ax[1].plot(M_tot, means, label = "mean log(p_hat(y|a))")
ax[1].plot(M_tot, logpdf(p_y_a(a0), y)*ones(size(N_tot)), "r-", label = "log(p(y|a))")
ax[1].set_xscale("log")
ax[1].grid(true)
ax[1].legend()
ax[1].set_xlabel("M")

ax[2].plot(M_tot, vars)
ax[2].set_xscale("log")
ax[2].set_yscale("log")
ax[2].grid(true)
ax[2].set_xlabel("M")
ax[2].set_ylabel("var log(p_hat(y|a))")

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.