Expected value of last Gamma RV in a sum I've got a sum of $X_i \sim \text{Gamma}(k, \theta)$ i.i.d. random variables. I'm trying to find the expected value of the final $X_i$ that takes the sum above a certain value, i.e., to find the value of 
$E[X_n | X_1 + X_2 + \cdots + X_{n-1} < 1000, X_1 + X_2 + \cdots + X_n > 1000]$ for integer $n > 1$.
 A: For simplicity, denote $X_n$ by $\xi$ and $X_1 + \cdots + X_{n - 1}$ by $\eta$. By assumption and the summation property of Gamma distribution, we know that $\xi \sim \Gamma(k, \theta)$ and $\eta \sim \Gamma((n - 1)k, \theta)$. In addition, $\xi$ and $\eta$ are mutually independent. Thus the problem reduces to find the conditional expectation $E[\xi | \eta < 1000, \xi + \eta > 1000]$. Note that unlike the common case, this "conditional expectation" is conditioning on an event instead of a random variable ($\sigma$-algebra), thus the result is a scalar instead of a random variable.
First, by the expectation formula for non-negative random variables, 
\begin{align}
& E[\xi | \eta < 1000, \eta + \xi > 1000] \\
= & \int_0^\infty P[\xi > t | \eta < 1000, \eta + \xi > 1000] dt \\
= & \frac{1}{P[\eta < 1000, \eta + \xi > 1000]}\int_0^\infty P[\xi > t, \eta < 1000, \eta + \xi > 1000] dt \tag{1}
\end{align}

Edit: Just realized an analytical solution is possible by means of Gamma function and Beta function when the shape parameter $k$ is a positive integer. To save some typing, let's assume without generality that the scale parameter $\theta \equiv 1$, and denote $1000$ by $\alpha$. Moreover, let 
\begin{align}
f_\xi(x) & = \frac{1}{\Gamma(k)}e^{-x}x^{k - 1} \equiv c_1 e^{-x}x^{k - 1} \quad (x > 0), \\
f_\eta(y) & = \frac{1}{\Gamma((n - 1)k)}e^{-y}y^{(n - 1)k - 1} \equiv c_2 
e^{-y}y^{m - 1} \quad (y > 0) \\
\end{align}
denote density functions of $\xi$ and $\eta$. Then 
$$P[\eta < \alpha, \xi + \eta > \alpha] = c_1c_2\int_0^\alpha\left[\int_{\alpha - y}^\infty e^{-x}x^{k - 1}dx\right] 
e^{-y}y^{m - 1}dy. \tag{2}$$
Make substitution $t = x - (\alpha - y)$ in the inner integral of $(2)$, it follows that 
\begin{align}
& \int_{\alpha - y}^\infty e^{-x}x^{k - 1} dx \\
= & e^{y - \alpha} \int_0^\infty e^{-t}[t + (\alpha - y)]^{k - 1} dt \\
= & e^{y - \alpha} \sum_{i = 0}^{k - 1} \binom{k - 1}{i}(\alpha - y)^i\int_0^\infty e^{-t}t^{(k - i) - 1}dt \\
= & e^{y - \alpha} \sum_{i = 0}^{k - 1} \binom{k - 1}{i}(\alpha - y)^i\Gamma(k - i).
\end{align}
Therefore, the right hand side of $(2)$ equals to 
\begin{align}
& c_1c_2 e^{-\alpha}\sum_{i = 0}^{k - 1}\binom{k - 1}{i}\Gamma(k - i)\int_0^\alpha
y^{m - 1}(\alpha - y)^i dy \\
= & c_1c_2e^{-\alpha} \sum_{i = 0}^{k - 1}\binom{k - 1}{i}\Gamma(k - i)\alpha^{m - i}\int_0^1 z^{m - 1}(1 - z)^{i + 1 - 1} dz \\
= & c_1c_2e^{-\alpha} \sum_{i = 0}^{k - 1}\binom{k - 1}{i}\Gamma(k - i)\alpha^{m - i}B(m, i + 1) \\
= & c_1c_2e^{-\alpha} \sum_{i = 0}^{k - 1}\binom{k - 1}{i}\Gamma(k - i)\alpha^{m - i}\frac{\Gamma(m)\Gamma(i + 1)}{\Gamma(m + i + 1)} \\
= & c_1c_2\alpha^me^{-\alpha}\sum_{i = 0}^{k - 1}\frac{(k - 1)!(m - 1)!}{(m + i)!\alpha^i} = \alpha^me^{-\alpha}\sum_{i = 0}^{k - 1}\frac{1}{(m + i)!\alpha^i}.
\end{align}
Now let's treat the numerator of $(1)$, careful calculation shows that 
\begin{align}
& \int_0^\infty P[\xi > t, \eta < \alpha, \xi + \eta > \alpha] dt \\
= & \int_0^\infty \int_t^\infty P_\eta[\alpha - x < \eta < \alpha]f_\xi(x)dx  dt \\
= & \int_0^\infty\int_t^\infty P[\eta < \alpha]f_\xi(x) dx dt - \int_0^\alpha\int_t^\alpha P[\eta < \alpha - x]f_\xi(x) dx dt \\
= & P[\eta < \alpha]E[\xi] - \int_0^\alpha\int_t^\alpha P[\eta < \alpha - x]f_\xi(x) dx dt \\
= & kP[\eta < \alpha] - \int_0^\alpha\int_t^\alpha P[\eta < \alpha - x]f_\xi(x) dx dt.
\end{align}
The latter term of the above integral can be expanded as:
\begin{align}
& c_1c_2\int_0^\alpha \int_t^\alpha\left[\int_0^{\alpha - x} e^{-y}y^{m - 1}dy\right]e^{-x}x^{k - 1}dxdt \\
= & c_1c_2\int_0^\alpha\left[\int_0^{\alpha - y}e^{-x}x^{k}dx\right]e^{-y}y^{m - 1}dy \\
= & c_1c_2\int_0^\alpha\left[\Gamma(k + 1) - \int_{\alpha - y}^\infty e^{-x}x^{k}dx\right]e^{-y}y^{m - 1}dy \\
= & c_1k!P[\eta < \alpha] - c_1c_2\int_0^\alpha\left[\int_{\alpha -y}^\infty e^{-x}x^{k}dx\right]e^{-y}y^{m - 1}dy \\
= & kP[\eta < \alpha] - c_1c_2\int_0^\alpha\left[\int_{\alpha -y}^\infty e^{-x}x^{k}dx\right]e^{-y}y^{m - 1}dy.
\end{align}
Therefore
$$\int_0^\infty P[\xi > t, \eta < \alpha, \xi + \eta > \alpha] dt =
c_1c_2\int_0^\alpha\left[\int_{\alpha -y}^\infty e^{-x}x^{k}dx\right]e^{-y}y^{m - 1}dy. \tag{3} $$
Fortunately, $(3)$ is almost the same as $(2)$, except for $k$ in $(2)$ will be replaced by $k + 1$ in $(3)$, whence a direct analogy gives 
$$\int_0^\infty P[\xi > t, \eta < \alpha, \xi + \eta > \alpha] dt =
\alpha^me^{-\alpha}\sum_{i = 0}^{k}\frac{1}{(m + i)!\alpha^i}. $$
Finally, we can conclude that 
$$E[\xi | \eta < 1000, \xi + \eta > 1000] = \frac{\sum_{i = 0}^{k}\frac{1}{((n - 1)k + i)!1000^i}}{\sum_{i = 0}^{k - 1}\frac{1}{((n - 1)k + i)!1000^i}}.$$
