Confidence Interval for Inverse Gamma Distribution

Question

I would like to understand if there exists any method to find confidence interval for the parameters of inverse gamma distribution.

Answer 1

Some examples with R (and various extension packages). First simulate some test data:

library(extraDistr)
library(MASS)
library(fitdistrplus)  

set.seed(7*11*13) # for reproducibility
x  <-  rinvgamma(200, 1, 1)

Then using functions from fitdistrplus package:

> xfit1 <-   fitdist(x, "invgamma", start=list(alpha=1.1,  beta=0.9))
xfit1
confint(xfit1)
> 
Fitting of the distribution ' invgamma ' by maximum likelihood 
Parameters:
      estimate Std. Error
alpha 1.075588 0.09531133
beta  1.089889 0.12185770
>  
         2.5 %   97.5 %
alpha 0.8887813 1.262395
beta  0.8510528 1.328726

But it does not seem like likelihood profiling is implemented for these objects, so this is only the usual asymptotic confidence interval, calculated from information in

 vcov(xfit1)
 coef(xfit1)
            alpha        beta
alpha 0.009084249 0.009205023
beta  0.009205023 0.014849299
>  
  alpha     beta 
1.075588 1.089889 

So for likelihood profiling we go for the package bbmle, which also is reasonably easy to use:

library(bbmle)
dat  <-  data.frame(x)


> xfit2  <-  mle2( x ~ dinvgamma(alpha=alpha, beta=beta), start=list(alpha=0.9, beta=0.9), data=dat, lower=c(alpha=0.1, beta=0.1), upper=c(alpha=10, beta=10), method="L-BFGS-B" )
> xfit2

Call:
mle2(minuslogl = x ~ dinvgamma(alpha = alpha, beta = beta), start = list(alpha = 0.9, 
    beta = 0.9), method = "L-BFGS-B", data = dat, lower = c(alpha = 0.1, 
    beta = 0.1), upper = c(alpha = 10, beta = 10))

Coefficients:
    alpha      beta 
0.7987156 0.8081280 

Log-likelihood: -248.01 
> confint(profile(xfit2))
          2.5 %   97.5 %
alpha 0.6239403 1.006797
beta  0.5711298 1.096452

Note that the two intervals (asymptotic and based on likelihood profiling) are rather different, and for these data the asymptotic ones looks better. These differences look too large to be a purely numeric issue, so warrants more investigation.

Answer 2

Since inverse gamma distributions are so often used in Bayesian Inference, another approximate finite sample inference approach is to use MCMC or Gibbs sampling to draw from a posterior using an uninformative prior to obtain credible intervals. Whilst credible != confidence, most agree the approach yields approximately equivalent inference when using non-informative priors.

Using the simulation suggestion from @KjetilBHalvorsen, I generate the same X and fit a BUGS model with:

model {
  for(i in 1:200) {
    invx[i] <- 1/x[i]
    invx[i] ~ dgamma(shape, scale)
  }
    shape ~ dgamma(0.1, 0.1)
    scale ~ dgamma(0.1, 0.1)

}

To obtain the following posterior distribution samples, which have 2.5 and 97.5 quantiles given by

Inference for Bugs model at "cat.txt", fit using WinBUGS,
 2 chains, each with 5000 iterations (first 2500 discarded), n.thin = 5
 n.sims = 1000 iterations saved
            mean    sd    2.5%     25%     50%     75%   97.5%  Rhat n.eff
shape      1.070 0.096   0.887   1.009   1.067   1.129   1.280 1.000  1000
scale      1.084 0.122   0.860   1.001   1.078   1.164   1.329 1.002   790
deviance 396.018 1.994 394.100 394.600 395.400 396.800 401.800 1.005   520

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 2.0 and DIC = 398.0
DIC is an estimate of expected predictive error (lower deviance is better).

Which agrees somewhat with the maximum likelihood approach used above.