I'm reading about Laplace's method for approximating the normalizing constant of a probability density function (which often is hard to compute explicitly) (MacKay, Information Theory, Inference, and Learning Algorithms, Chapter 27 – Laplace's Method).
I believe to have understood the concept. To visualize the individual steps involved I tried to approximate the normalizing constant of a peaked Gamma-distribution:
$$P(x)=\frac{1}{Z}(\frac{x}{s})^{c-1}exp(- \frac{x}{s})$$
where $Z=\Gamma(c)*s$,
and $s=1$, $c=3$.
Thus $Z=3!*1=6$.
The unnormalized distribution (w/o the $Z$) is:
$$P^*(x)=(\frac{x}{s})^{c-1}exp(- \frac{x}{s})$$
I am now assuming that I only know $P^*$ and that I'm trying to find an approximation for $Z$.
To do this I took the logarithm of $P^*$ and determined the location of its maximum $x_0$:
$$ln[P^*(x)]=(c-1)*ln(\frac{x}{s})-\frac{x}{s}$$
$$\frac{\partial ln[P^*(x)]}{\partial x} = (c-1)\frac{1}{x}-\frac{1}{s}=0 \implies x_0=2$$
Then I Taylor-expand up to the second derivative around the maximum and use it as an approximation to $P^*$:
$$ln[P^*(x)] \approx ln[P^*(x_0)]+\frac{1}{2}*d*(x-x_0)^2$$
$$d=\frac{\partial^2 ln[P^*(x_0)]}{\partial x^2} = -(c-1)\frac{1}{x^2}$$
When exponentiating I get an approximation $Q^*(x)$ to $P^*(x)$ which is of Gaussian form:
$$Q^*(x)=P^*(x_0)exp(\frac{d}{2}(x-x_0)^2)$$
which has normalizing constant:
$$Z_Q=\int Q^*(x)dx = P^*(x_0) \sqrt{\frac{2\pi}{-d}} $$
Yet, when evaluating $Z_Q\approx 1.91$, so threefold smaller than $Z$.
Can someone point out what I'm doing wrong?
Best Answer
The correct normalization constant is $Z = s \cdot \Gamma(c)=(c-1)! = 2$ so the Laplace method isn't too bad.