Solution (source code)

= Solution

Let $s=\sigma^2$ and $\vartheta=(s,\tau_1,\ldots,\tau_m)$. With <conditional independence> of the binary observations given the latent <function>, define
$$
L(\vartheta)=\int\ell(F)\varphi_{C_\vartheta}(F)\,dF,
\qquad \ell(F)=\prod_{i=1}^n\Phi((2y_i-1)F_i).
$$
The <prior distribution> is specified on the precision $s^{-1}$, not on $s$. Its <Jacobian determinant> gives
$$
p(\vartheta)=s^{-2}e^{-1/s}\prod_{j=1}^m e^{-\tau_j},\qquad s>0,\ \tau_j>0,
$$
and the marginal <posterior density> is proportional to $p(\vartheta)L(\vartheta)$. Working instead with precision or logarithmic coordinates is possible, provided all corresponding <Jacobian determinants> are included.

Choose an <importance sampling> density $g_\vartheta$ covering the support of $\ell(F)\varphi_{C_\vartheta}(F)$, draw $F_1,\ldots,F_M$ independently from it, and use the nonnegative <unbiased likelihood estimator>
$$
\widehat L(\vartheta)=\frac1M\sum_{j=1}^M
\frac{\ell(F_j)\varphi_{C_\vartheta}(F_j)}{g_\vartheta(F_j)}.
$$
Each summand integrates to $L(\vartheta)$, proving unbiasedness. A simple always-valid choice is the latent Gaussian <prior distribution>, $g_\vartheta=\varphi_{C_\vartheta}$, when the <covariance matrix> is nonsingular. Then the estimator is just the average of the bounded values $\ell(F_j)\in(0,1)$. For duplicate inputs, sample only the distinct latent values and multiply all likelihood contributions at each common input. A Gaussian approximation to the latent <posterior distribution>, or its mixture with the Gaussian <prior distribution>, can reduce relative <variance> while keeping complete support. Its approximation does not replace the exact weights.

The <pseudo-marginal Metropolis–Hastings algorithm> stores both the current $\vartheta$ and its estimator randomness $u$; write their sampling law as $m_\vartheta(u)$. Propose $\vartheta'\sim q(\vartheta,\cdot)$, independently generate $u'\sim m_{\vartheta'}$, compute $\widehat L(\vartheta',u')$, and accept the pair with
$$
\boxed{\alpha=1\wedge
\frac{p(\vartheta')\widehat L(\vartheta',u')q(\vartheta',\vartheta)}
{p(\vartheta)\widehat L(\vartheta,u)q(\vartheta,\vartheta')}.}
$$
Here $q(\vartheta,\vartheta')$ denotes the proposal <probability density function> from $\vartheta$ to $\vartheta'$. On rejection retain both the old parameter and the old likelihood estimate. Initialize with a positive estimate. One must not independently refresh the denominator after rejection.

Correctness follows from the extended <posterior density>
$$
\widetilde\pi(\vartheta,u)\propto
p(\vartheta)\widehat L(\vartheta,u)m_\vartheta(u).
$$
The pair proposal is $q(\vartheta,\vartheta')m_{\vartheta'}(u')$. Its <Metropolis–Hastings acceptance probability> simplifies to the displayed ratio, since the auxiliary densities cancel. Hence <detailed balance> holds for $\widetilde\pi$. Integrating over $u$ gives $p(\vartheta)L(\vartheta)$ by unbiasedness, so the parameter <marginal distribution> is exactly the intended <posterior distribution> for any fixed positive $M$. This is exactness of the stationary target, not independent exact samples or a guarantee of rapid convergence. Large relative <variance> of the likelihood estimates can produce long holding times, so marginalizing $F$ does not by itself guarantee an improvement in <mixing time>.