Previously we had the following hypothesis test for H0:a=a0 vs Ha:a≠a0: reject the null hypothesis if
so the acceptance region of the test is
to get a confidence interval we need to "invert the test". This means solving the double-inequality in the acceptance region for a (which now replaces a0). The next graph shows the left and the right side of the double-inequality as a function of a, (blue=α/2, red=1-α/2, n=50, 
=0.4)
and we get an interval (0.48, 1.45)
Formally we need to solve the equation beta2.ex(1) does the calculations.
How about the mle? This one is easy: Bayesian Interval Which Interval is Better?
First of frequentist confidence intervals and Bayesian credible intervals can not really be compared directly. We will just compare the two confidence intervals.
First we need to check that the two method yield proper confidence intervals, that is that they have coverage. So if we calculate a 90% interval it really contains the true parameter 90% of the time. This is true for the LRT interval because we could find the exact distribution of the LRT statistic and invert the interval analytically. The other test is a large sample test and needs to be checked via simulation. here are some results (red=MM, blue=LRT) How can we choose between the two? We need a measure of performance for a confidence interval. Here is a possible choice: the mean length of the interval E[U(X)-L(X)]. In the case of the MM method this will need to be found via simulation. As for the LRT method, we have Above we used the fact that aT has distribution that does not depend on the parameter a. We then found the interval by dividing the error probability α equally on the left and the right. But is that the best option? If we are interested in shortest-length intervals, can we find the intervals that are optimal? That is can we find L(T) and U(T) such that
minimize E[U(T)-L(T)] subject to P(L(T)<a<U(T))=1-α
Let's set q(z)=qgamma(z,n,1) and let's consider intervals of the form L(T)=q(z1)/T and U(T)=q(1-z2)/T. Clearly we need z1+z2=α , so we have z2=1-α+z1 or we just write
L(T)=q(z)/T and U(T)=q(1-α+z)/T
This has length How much better are these intervals? Here are the mean lengths (for a=1)
for z=zα/2 and z=z1-α/2. Let x=, then
so now we have a cubic equation, which we can solve. In R this is done by the routine cuberoots. Of course there are generally three roots, but one of them is negative and the other two are our solutions. Note that zα/22=qnorm(1-α/2)2=qnorm(α/2)2=z1-α/22 , so we only need to solve one equation.
beta2.ex(2) does the work.
Previously we saw that if we use a prior Exp(1) we get a posterior distribution a|x~Γ(n+1,1/(T+1)) using this we can find the equal-tail probability interval by solving
which is (almost) the same as the confidence interval based on the LRT test.
so we see that the MM method has some over-coverage for small n. Over-coverage is acceptable but not really desirable.
Here is a graph of the mean lengths:
This can not be done analytically, and so again we need to resort to a numerical solution. In beta2.ex(3) i use a simple "grid-search" algorithm. It turns out that assigning a smaller part of α to the left side leads to shorter intervals.
For example if n=2 (and a=1) the mean length of the equal-tail intervals is 5.33 but for the "optimal" length intervals it is 4.72, an improvement of about 10%