5

I'm reading Example 3.16 of Robert & Casella's Monte Carlo Statistical Methods. It uses a Laplace approximation for approximating an integral related with the Gamma distribution namely $$\int_a^b\frac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}e^{-\frac{x}{\beta}}dx$$

I'm wondering where $\Gamma(\alpha)\beta^\alpha$ does go in the end result. Also, the book didn't mention too much about importance of $n$, can anyone give some explanation on this regard?

enter image description here

Xi'an
  • 90,397
  • 9
  • 157
  • 575
ZHU
  • 505
  • 3
  • 11
  • 1
    You're going to need to quote more of the book. Looking only at your question here (I don't have the book), I have no idea where it went in the end result. I didn't know it went anywhere until I read your 2nd paragraph, eg. – gung - Reinstate Monica Jun 29 '18 at 14:33
  • 1
    @gung I've added the page – ZHU Jul 01 '18 at 07:26
  • @gung can you re-open this question? – ZHU Jul 01 '18 at 07:27
  • @Xi'an thanks for your reply! I was wondering if you can talk about why $$n$$ needs to be included in general? It seems to me it didn't affect the Taylor expansion and all that. – ZHU Jul 25 '18 at 06:45
  • The Taylor expansion is only valid in the neighbourhood of the value at which the expansion takes place. The farther one goes the less appropriate it gets. When considering the integral example, the only infinitesimal one can spot in the quantity is the length of the integral. When $b-a$ goes to zero, the expansion gets exact. – Xi'an Jul 25 '18 at 20:55

1 Answers1

8

This example is indeed rather poorly conducted and full of typos, apologies from one author!!!

First, there is no genuine $n$ (or related sample size) in the picture so the Laplace approximation cannot be universally adequate. The only thing that can replace $n$ in the integral $$ \int_a^b \dfrac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}e^{-x/\beta}\text{d}x$$ is $(b-a)^{-1}$, in the sense that the approximation becomes better and better as $(b-a)$ goes to zero. This is under the provision that the mode $\hat{x}_\theta$ belongs to the interval $(a,b)$.

Second, as you correctly wondered, the constant $\Gamma(\alpha)\beta^\alpha$ in the denominator has truly and inexplicably (!) been forgotten in the final formula, which should thus be \begin{align}\int_a^b &\dfrac{x^{\alpha-1}}{\Gamma(\alpha)\beta^\alpha}e^{-x/\beta}\text{d}x\,\approx \\ &\dfrac{\hat{x}_\theta^{\alpha-1}e^{-\hat{x}_\theta/\beta}}{\Gamma(\alpha)\beta^\alpha}\sqrt{\frac{2\pi \hat{x}_\theta^2}{\alpha-1}}\left\{\Phi\left(\sqrt{\frac{\alpha-1}{2\pi \hat{x}_\theta^2}} (b-\hat{x}_\theta)\right)-\Phi\left(\sqrt{\frac{\alpha-1}{2\pi \hat{x}_\theta^2}} (a-\hat{x}_\theta)\right)\right\}\end{align}

Third, the formula at the top of page 110 is missing two minus ($-$) signs, as it should in truth be $$h(x)\approx -\frac{\hat{x}_\theta}{\beta}+(\alpha-1)\log \hat{x}_\theta-\frac{\alpha-1}{2\hat{x}_\theta^2}(x-\hat{x}_\theta)^2$$

With these multiple corrections, the table comparing the Laplace version with the exact coverage can be reconstructed, for instance via the following R code:

xop <- function(al,be){(al-1)*be}

lapl <- function(al,be,a,b){
 hatx=xop(al,be)
 hatx^{al-1}*exp(-hatx/be)*sqrt(2*pi)*hatx/sqrt(al-1)*
 (pnorm(sqrt(al-1)*(b-hatx)/hatx)-pnorm(sqrt(al-1)*(a-hatx)/hatx))/
 be^al/factorial(al-1)
 }

xact <- function(al,be,a,b){
  pgamma(b,al,1/be)-pgamma(a,al,1/be)}

> xact(5,2,7,9)
[1] 0.1933414
> lapl(5,2,7,9)
[1] 0.1933507

Once again, my most sincere apologies for this poor attention to details in this example!

Xi'an
  • 90,397
  • 9
  • 157
  • 575
  • 2
    (+1) Answers don't get much more authoritative than this! – Sycorax Jul 18 '18 at 18:03
  • 2
    It's pretty great that you can post a question on CV.SE and the co-author of the statistic book you are asking about writes the answer! – Ben Jul 22 '18 at 10:41