Argthtjtr

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
  

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.

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=logfrac{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*}