Simulating the Posterior Density of a Transformed Parameters

Question

I am reviewing an example (p. 180-181, Example 11.3 and 11.4) from All of Statistics by Larry Wasserman. The example intends to illustrate that the posterior can be found analytically and can be approximated by simulation approach:

Let $X_1,...X_n \sim Bernoulli(p)$ and $f(p) = 1$ so that $p|X_1,...,X_n \sim Beta(s+1, n-s+1)$ with $s = \sum_{i=1}^{n} x_i $. Let $\psi = log(p/(1-p))$. Find the posterior density of $\Psi|X_1,...,X_n$.

By analytical approach, it is shown that the exact pdf of the posterior is as follows:

$h(\psi|x_1,...,x_n) = H'(\psi|x_1,...,x_n) = \frac{\Gamma(n+2)}{\Gamma(s+1)\Gamma(n-s+1)}(\frac{e^\psi}{1+e^\psi})^{s}(\frac{1}{1+e^\psi})^{n-s+2}$

To approximate the above exact pdf, one can perform simulation approach as follows:

1. Draw $P_1,...,P_B \sim Beta(s+1, n-s+1)$
2. Let $\psi_i = log(P_i / (1 - P_i))$ for $ i = 1,...,B.$ Now $\psi_1,...,\psi_B$ are IID draws from $h(\psi|x_1,...,x_n).$
3. Plot the histogram for $\psi.$

Based on the above procedures, I tried to do a simulation experiment to see if the simulated result is actually closed to the exact pdf. The plots below are my result:

The upper graph is the exact pdf derived by analytical approach while the lower graph is the simulated histogram. Their shapes are pretty similar, but they are different in the heights. Particularly, the peak of the simulated histogram is greater than 1 even after the normalisation, while the peak of the exact pdf is below 1. I am trying to figure out the reason for such discrepancy. Could anyone provide me any clues? Thanks!

• You can find my whole work in the link below: Full Simulation Result

Answer 1

To anyone coming from Google: run the code provided by the asker prior to running my code.

The PDF you provide does not seem to integrate to 1 (it integrates to 33/68 for the seed you provide, as confirmed using Mathematica).

We can numerically integrate the PDF in R and see if it looks right. Unfortunately, numerical stability complicates this for the PDF as you have it written:

> integrate(psi_pdf, -Inf, Inf)
Error in integrate(psi_pdf, -Inf, Inf) : non-finite function value

This is due to the combinatorial gamma functions and exponentials present in the function leading to numerical issues. This may be remedied by rewriting the pdf in logspace, then exponentiating at the end:

> stable_pdf <- function(psi) {
>+ lgamma(n+2) - lgamma(s+1) - lgamma(n-s+1) + s * (psi - log(1+exp(psi))) + (n - s + 2) * (-log(1+exp(psi)))
>+ }

We see that this matches up with the old pdf if logged:

> stable_pdf(0)
[1] -5.860918
> log(psi_pdf(0))
[1] -5.860918

And can see that the integral is far from 1:

> integrate(function(psi) exp(stable_pdf(psi)), -Inf, Inf)
0.4852941 with absolute error < 1.5e-05

Just as your graphics suggest, the PDF needs to be multiplied by approximately 2 to have total measure 1. Without access to the book example, I cannot say where the error occurs in its derivation.

Answer 2

The density of $$p|X_1,...,X_n \sim \text{Beta}(s+1, n-s+1)$$ is $$f(p|X_1,...,X_n) = \frac{\Gamma(n+2)}{\Gamma(s+1)\Gamma(n-s+1)} p^s (1-p)^{n-s}$$ Hence the density of $$\psi=\log\frac{p}{1-p}$$ is $$h(\psi|X_1,...,X_n) = f(1/1+e^{-\psi})\overbrace{\left|\frac{\text{d}p}{\text{d}\psi}\right|}^\text{Jacobian}$$that is \begin{align*}h(\psi|X_1,...,X_n) &= \frac{\Gamma(n+2)}{\Gamma(s+1)\Gamma(n-s+1)} \frac{e^{\psi s}}{(1+e^\psi)^n} \frac{e^{\psi}}{(1+e^\psi)^2}\\ &=\frac{\Gamma(n+2)}{\Gamma(s+1)\Gamma(n-s+1)} \frac{e^{\psi (s+1)}}{(1+e^\psi)^{n+2}}\end{align*}