The Harmonic Mean of the Likelihood, n sample normal-normal case

Question

This is very related to Radford Neal's awesome blogpost, but his example uses a single data point x, and I want to know how to replicate the results using $n$ data points. This is essentially what my question will "boil down" to. Below I'll try to describe what I've done, but its not essential to read if you already know what I'm asking for.

Just as he explains, we have a normal-normal case, so we can compute the posterior distribution analytically. If $X \sim N(\theta, 1)$ with $\theta \sim N(0,1)$, then this time (since $n>1$) we get that the

$$p(\theta|x) \propto \exp\left\{ -\frac{1}{2}\left(\frac{1+1/n}{1}\right)\left(\theta-\frac{\bar{x}}{1+1/n}\right)\right\}, $$

i.e. the posterior distribution of $\theta$ is $$ N\left( \frac{\bar{x}}{1+1/n} \, , \,\frac{1/n}{1+1/n} \right).$$

My procedure should look something like:

1. draw k sample points from the posterior distribution of theta
2. for each point, compute the likelihood.
3. compute reciprocal for each of the likelihoods
4. compute the mean of them
5. compute the reciprocal of the mean, and let this value be the estimator for the marginal likelihood p(y)

Step 1 is easy, say that $n=50$, we can then sample $k$ points directly from the posterior distribution using

    data <- rnorm(50)
    n <- length(data)
    post_mean <- mean(data)/(1+1/n) 
    post_sd   <- sqrt((1/n) / (1+1/n))
    post_sample <- rnorm(k, post_mean, post_sd)

Step 3-5 are straight forward. I'm stuck at stage 2.

Now this is what the function from the article above looks like.

    harmonic.mean.marg.lik <- function (x, s0, s1, n)
    { post.prec <- 1/s0^2 + 1/s1^2
    t <- rnorm (n, (x/s1^2)/post.prec, sqrt(1/post.prec))
    lik <- dnorm(x,t,s1)
    1/mean(1/lik)
    }

Clearly, it is lik <- dnorm(x,t,s1) which bothers me. His t is my post_sample, his s0 is $1$, and in my case his s1 is $1/n$. His x however is just a single data point. And I think maybe I'm confused about what R is doing if I were to throw in a vector with data. Because I do get a result, but I am not sure what is going on. For example, if I set $k=10000$, then this would be my attempt

    > k = 10000
    > data <- rnorm(1000)
    > n <- length(data)
    > post_mean <- mean(data)/(1+1/n) 
    > post_sd   <- sqrt((1/n) / (1+1/n))
    > post_sample <- rnorm(k, post_mean, post_sd)
    > lik <- dnorm(data, post_sample, 1/n )
    > 1/mean(1/lik)
    [1] 0

But I don't know if what I'm doing is right or not.

Answer 1

His x however is just a single data point.

Thanks to the sufficiency of $\bar{x}$ for a N$(\theta,1)$ distribution, the likelihood of the sample is (proportional to) the normal density at $\bar{x}$ for a $1/n$ variance. There is therefore no need to compute the likelihood as a product of normal densities at each point $x_i$ of the sample.

When calling

lik <- dnorm(data, post_sample, 1/n)

for data and post_sample of the same size, the returned vector is made of the $\varphi(\sqrt{n}(x_i-\theta_i)$'s which is not what you want. E.g.,

> dnorm(1:10,1:10,1/10)
 [1] 3.989423 3.989423 3.989423 3.989423 3.989423 3.989423 3.989423 3.989423
 [9] 3.989423 3.989423