Conditional expectation of a Gamma Distribution

Question

$X$ a Gamma Distributed random variable. Calculate $E[X\mid X \in [a,b]]$, where $a>0$, $b>0$. Is there a closed form solution for this, and if so, how can I calculate it?

Answer 1

What you are looking for is the expectation of a truncated gamma distribution. Formulas (11) and (13) in "A Right and Left Truncated Gamma Distribution with Application to the Stars" by Zaninetti (there is a pdf here) give you the formula you are looking for. Let $b$ denote the scale and $c$ the shape. Then

$$ E(X|X\in[x_\ell,x_u]) = b^2k\bigg(\Gamma\Big(1+c, \frac{x_\ell}{b}\Big)-\Gamma\Big(1+c,\frac{x_u}{b}\Big)\bigg), $$

where $\Gamma$ denotes the upper incomplete gamma function and

$$ k= \frac{c}{b\Gamma(1+c,\frac{x_\ell}{b})-b\Gamma(1+c,\frac{x_u}{b})+ e^{-\frac{x_u}{b}}b^{-c+1}x_u^c-e^{-\frac{x_\ell}{b}}b^{-c+1}x_\ell^c}.$$

I like verifying calculations like this using an R script, like this (note that pracma::incgam() switches the order of the two parameters of the upper incomplete gamma function compared to the formulation I took from the paper):

require(pracma)

shape <- 2
scale <- 3
aa <- 1
bb <- 4

set.seed(1)
foo <- rgamma(1e6,shape,scale=scale)
mean(foo[foo>aa & foo<bb])

kk <- shape/
    (scale*incgam(aa/scale,1+shape)-scale*incgam(bb/scale,1+shape)+
        exp(-bb/scale)*scale^(-shape+1)*bb^shape-exp(-aa/scale)*scale^(-shape+1)*aa^shape)
scale^2*kk*(incgam(aa/scale,1+shape)-incgam(bb/scale,1+shape))

The call to mean() and the last command give the same result up to noise, also for other values of the parameters.

Answer 2

The formula for the expected value of a gamma random variable (with shape parameter $\alpha$ and scale parameter $\beta$) constrained to an interval $\left[ a,b \right]$ can be expressed as $$ E \left[ X \ | \ a<X<b \ \right]=\frac{\alpha \beta \left[ P \left( \alpha+1,\frac{b}{\beta} \right) - P \left( \alpha+1,\frac{a}{\beta} \right) \right] }{P \left( \alpha,\frac{b}{\beta} \right) - P \left( \alpha,\frac{a}{\beta} \right)} \ , $$ where the function $P \left( \alpha,x \right)$ is the lower incomplete gamma function defined by $$P \left( \alpha,x \right) = {\frac{1}{\Gamma \left( \alpha \right)}} \int_0^x t^{\alpha-1} \ e^{-t} \ dt$$

Here is R code to illustrate:

# Assign shape and scale
alpha <- 3
beta <- 2

# Simulate
x <- rgamma(5000000,shape=alpha,scale=beta)

# Constrain to interval
a <- 3
b <- 5
y <- x[x > a & x < b]

# Get sample mean
mean(y)

# Calculate expected value
numer <- alpha*beta*(pgamma(b/beta,shape=alpha+1,scale=1)-pgamma(a/beta,shape=alpha+1,scale=1))

denom <- pgamma(b/beta,shape=alpha,scale=1)-pgamma(a/beta,shape=alpha,scale=1)

expec <- numer/denom

expec

Here's the output:

> mean(y)
[1] 4.002289


> expec
[1] 4.002089