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

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$$?

• Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
– Community Bot
Jan 9 at 16:48

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.

• Yours is indeed a brilliant solution to a very general problem. When working with asymmetric density functions, the mean and variance alone cannot help find the desired bounds provided a given total probability. As such, your answer applies to numerous situations in Bayesian workflows. Jan 10 at 10:33
• Thanks! If you think this answers the question, please accept it, as it will let future users of the site who may have similar questions that there is a helpful answer here. Jan 10 at 16:00