# Bias in unbiased pseudo-marginal estimation

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$$, The distribution of the estimates is shown below

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))")