Eliciting a Gamma informative prior in a Gamma–Poisson Bayesian problem

Question

1. I employ the Gamma–Poisson conjugate family for my statistical model.
2. I want to use an informative prior.
3. From theory, I know that the values of the Gamma-distributed random variable lie within the interval $[x_0,x_1]$.
4. Using the shape-and-scale characterization, I need to find the values of $k$ and $\theta$ such that the corresponding Gamma density contains a ninety-five percent probability between $x_0$ and $x_1$.
5. Using Gamma's CDF, I get the two conditions: $$\frac{1}{\Gamma(k)} \gamma \left( k, \frac{x_0}{\theta}\right)=0.025$$ and $$\frac{1}{\Gamma(k)} \gamma \left( k, \frac{x_1}{\theta}\right)=0.975$$ that need to be solved for $k$ and $\theta$. Both conditions involve $\gamma$, the lower incomplete Gamma function.

How can I solve for $k$ and $\theta$?

Answer 1

What you have is a two-dimensional nonlinear root finding problem - two-dimensional because you have two variables ($k$ and $\theta$) and two equations. It's a nice problem because both the functions are continuous and differentiable in both the variables, and I'm pretty sure there's only one root.

One way to solve well-behaved problems like this, if you don't have any multidimensional root finding routines available to you, is to convert it to a minimization problem with an objective function such as:

$$ \min_{k,\theta} f(k,\theta) = (F(x_1;k,\theta)-0.025)^2 + (F(x_2;k,\theta)-0.975)^2$$

Multidimensional minimization algorithms abound; here's an implementation of this approach in R.

foo <- function(theta, x1, x2) {
  e1 <- pgamma(x1, theta[1], theta[2]) - 0.025
  e2 <- pgamma(x2, theta[1], theta[2]) - 0.975
  e1*e1 + e2*e2
}

true_parms <- c(2.5, 2.5)  # k, theta

x1 <- qgamma(0.025, true_parms[1], true_parms[2])
x2 <- qgamma(0.975, true_parms[1], true_parms[2])

c(x1, x2)  # the true lower and upper 2.5% and 97.5% values

[1] 0.1662423 2.5665004

optim(c(1,1), foo, x1=x1, x2=x2)  # Starting values are (1,1)
$par
[1] 2.500241 2.500304  # The results of the optimization routine

$value
[1] 8.802968e-11  # ... and other, marginally relevant, stuff.