Issues with the Bayes estimator of a conjugate normal hierarchy model

Someone asks for an instability issue in R's integrate program:

Hello everyone,

I am supposed to calculate the Bayes estimator of a conjugate normal hierarchy model. However, the Bayes estimator does not have a closed form,

The book "Theory of Point Estimation" claims that the numerical evaluation of  this estimator is simple. But my two attempts below both failed.

1. I tried directly using the integration routine in R on the numerator and denominator separately. Maybe because of the infinite domain, occasionally the results are far from reasonable.

2. I tried two ways of change of variables so that the resulting domain can be finite. I let


But the estimator results are very similar to the direct integration on the original integrand. More often than it should occur, we obtain quite large evaluation of the Bayes estimator, up to 10^6 magnitude.

I wonder if there is any other numerical integration trick which can lead to a more accurate evaluation.

I appreciate any suggestion.


-------------------------------------------
xxx
Some State University
-------------------------------------------

Well, what happens here? Her program have a part that says:

[Bayes(nu,p,sigma,xbar) is the ratio of both integrals, "f", and "g" are the integrals, f is the numerator, g the denominator, so Bayes = f/g]

Now, executing Bayes(2,10,1,9.3) fails:

> Bayes(2,10,1,9.3)
[1] 1477.394

, which is much greater than the expected approx. 8.


I tried this with the same program, integrate, to do this simple case (dnorm is the normal distribution density):

> integrate(dnorm,0,1)
0.3413447 with absolute error < 3.8e-15
> integrate(dnorm,0,10)
0.5 with absolute error < 3.7e-05
> integrate(dnorm,0,100)
0.5 with absolute error < 1.6e-07
> integrate(dnorm,0,1000)
0.5 with absolute error < 4.4e-06
> integrate(dnorm,0,10000000000)
0 with absolute error < 0



As we can see, the last try, with a very large value, fails miserably, value is 0 (instead of 0.5) and absolute error is negative.

The program "integrate" uses code that is part of supposedly "a Subroutine Package for Automatic Integration", as it is advertised, but it cannot anticipate everything -- and we hit an instability we cannot solve.

My suggestion was to use integrate(f,0,1) and integrate(g,0,1) always until we get results outside what is reasonable. In those cases, we should try integrate(f,0,.999) and integrate(g,0,.999) with as many nines as we can (I got problems with just .9999, that's why I wrote .999 there).

Of course, you can always try a different method. Since this function is well-behaved, any simple method could be good enough.