# Calculate mean and variance of conditional sum of Gammas

Find $Mean(Z)$ and $Var(Z)$, where:

$Z = \begin{cases}W + X_1 & \text{if } W \leq c \\ W + X_2 & \text{if } W > c \end{cases}$

where $\quad W \sim \Gamma(k_0, \theta_0), \quad X_1 \sim \Gamma(k_1, \theta_1), \quad X_2 \sim \Gamma(k_2, \theta_2)$ are independent, and $c$ is a constant $> 0$.

My code in mathematica using mathStatica is as follows:

f=((w^(Subscript[k, 0]-1)*Exp[-w/Subscript[\[Theta], 0]])/(Gamma[Subscript[k, 0]]*Subscript[\[Theta], 0]^(Subscript[k, 0])))*((Subscript[x, 1]^(Subscript[k, 1]-1)*Exp[-Subscript[x, 1]/Subscript[\[Theta], 1]])/(Gamma[Subscript[k, 1]]*Subscript[\[Theta], 1]^(Subscript[k, 1])))*((Subscript[x, 2]^(Subscript[k, 2]-1)*Exp[-Subscript[x, 2]/Subscript[\[Theta], 2]])/(Gamma[Subscript[k, 2]]*Subscript[\[Theta], 2]^(Subscript[k, 2])));
domain[f]={{w,-Infinity,Infinity},{Subscript[x,1],-Infinity,Infinity} ,{Subscript[x,2],-Infinity,Infinity}} && {Subscript[k, 0]>0,Subscript[k, 1]>0,Subscript[k, 2]>0,Subscript[\[Theta], 0]>0,Subscript[\[Theta], 1]>0,Subscript[\[Theta], 2] >0,g >0};
z=\[Piecewise]  w+Subscript[x,1] w<=c  w+Subscript[x,2] w>c;
mean=Expect[z,f]
var=Var[z,f]


However, mathematica (Mathstatica) is unable to calculate the mean and variance of this expression. Is there a better way to calculate the mean and variance of Z?

The following Monte carlo simulation shows:

Ddata=RandomReal[GammaDistribution[121,0.09],100000];


mean=11.91 and variance 1.10

• Start by removing Subscript from your variable names, which leads to all sort of troubles and makes your code unreadable. Then use GammaDistribution Distributed, Probability and TransformedDistribution Oct 6, 2015 at 13:34

You can define distributions and calculate probabilities, for instance of the condition $w < c$.

Probability[
w < c
]

Piecewise[{{1 - Gamma[k0, c/th0]/Gamma[k0], c > 0}}, 0]


Once distributions are defined you can easily calculate Mean and Variance.

Mean[
TransformedDistribution[
w + x1,
{
}
]
]

k0 th0 + k1 th1

Variance[
TransformedDistribution[
w + x1,
{
}
]
]

k0 th0^2 + k1 th1^2


Your distribution is conditional so you need to incorporate that condition using Conditioned

Mean@TransformedDistribution[
Conditioned[w + x1, w > c],
{
}
]

Piecewise[{{k0*th0 + k1*th1, c <= 0}}, k1*th1 + (th0*Gamma[1 + k0, c/th0])/Gamma[k0, c/th0]]


Or with Piecewise

Mean@TransformedDistribution[
Piecewise[{{w + x1, w > c}, {w + x2, w < c}}],
{
}]

Piecewise[{{k0*th0 + k1*th1, c <= 0}}, ((k0*th0 + k2*th2)*Gamma[k0] + (k1*th1 - k2*th2)*Gamma[k0, c/th0])/Gamma[k0]]


Assuming[And @@ Table[i > 0, {i, {c, k0, k1, k2, th0, th1, th2}}],
FullSimplify[
Variance@TransformedDistribution[
Piecewise[{{w + x1, w > c}, {w + x2, w < c}}],
{
}]
]]

(1/(Gamma[k0]^2))((k2 th2^2 +
k0 th0 (th0 + 2 k1 th1 - 2 k2 th2)) Gamma[
k0]^2 + (k1 th1 - k2 th2) Gamma[k0, c/
th0] (-2 th0 Gamma[1 + k0] + (-k1 th1 + k2 th2) Gamma[k0, c/
th0]) + Gamma[
k0] (2 th0 (-k1 th1 + k2 th2) Gamma[
1 + k0] + (k1 (1 + k1) th1^2 -
2 k1 k2 th1 th2 + (-1 + k2) k2 th2^2) Gamma[k0, c/th0] +
2 th0 (k1 th1 - k2 th2) Gamma[1 + k0, c/th0]))

N[
((k0 th0 + k2 th2) Gamma[k0] + (k1 th1 - k2 th2) Gamma[k0, c/th0])/
Gamma[k0] /.
{
k0 -> 121
, th0 -> 0.09
, k1 -> 100
, th1 -> 0.01
, k2 -> 5.76
, th2 -> 0.2083
, c -> 12
}]

12.0634

(k2 th2^2 Gamma[
k0]^2 - (th0 Gamma[1 + k0] + (k1 th1 - k2 th2) Gamma[k0, c/
th0])^2 +
Gamma[k0] (th0^2 Gamma[
2 + k0] + (k1 (1 + k1) th1^2 -
2 k1 k2 th1 th2 + (-1 + k2) k2 th2^2) Gamma[k0, c/th0] +
2 th0 (k1 th1 - k2 th2) Gamma[1 + k0, c/th0]))/Gamma[k0]^2
/. {
k0 -> 121
, th0 -> 0.09
, k1 -> 100
, th1 -> 0.01
, k2 -> 5.76
, th2 -> 0.2083
, c -> 12
}]

1.1153203656002