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

  • $\begingroup$ 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. $\endgroup$
    – Community Bot
    Jan 9 at 16:48

1 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)
[1] 2.500241 2.500304  # The results of the optimization routine

[1] 8.802968e-11  # ... and other, marginally relevant, stuff.
  • $\begingroup$ 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. $\endgroup$
    – Valerio
    Jan 10 at 10:33
  • 1
    $\begingroup$ 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. $\endgroup$
    – jbowman
    Jan 10 at 16:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.