1
$\begingroup$

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];
T1data=RandomReal[GammaDistribution[100,0.01],100000];
T2data=RandomReal[GammaDistribution[5.76,0.2083],100000];
Adata=MapThread[If[#1<12,#1+#2,#1+#3]&,{Ddata,T1data,T2data}];
N[Mean[Adata]]
N[Sqrt[Variance[Adata]]]

mean=11.91 and variance 1.10

$\endgroup$
  • 1
    $\begingroup$ 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 $\endgroup$ – rhermans Oct 6 '15 at 13:34
3
$\begingroup$

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

Probability[
 w < c
 , Distributed[w, GammaDistribution[k0, th0]]
 ]
Piecewise[{{1 - Gamma[k0, c/th0]/Gamma[k0], c > 0}}, 0]

Mathematica graphics

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

Mean[
 TransformedDistribution[
  w + x1,
  {
   Distributed[w, GammaDistribution[k0, th0]]
   , Distributed[x1, GammaDistribution[k1, th1]]
   }
  ]
 ]
k0 th0 + k1 th1
Variance[
 TransformedDistribution[
  w + x1,
  {
   Distributed[w, GammaDistribution[k0, th0]]
   , Distributed[x1, GammaDistribution[k1, th1]]
   }
  ]
 ]
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],
  {
   Distributed[w, GammaDistribution[k0, th0]]
   , Distributed[x1, GammaDistribution[k1, th1]]
   }
  ]
Piecewise[{{k0*th0 + k1*th1, c <= 0}}, k1*th1 + (th0*Gamma[1 + k0, c/th0])/Gamma[k0, c/th0]]

Mathematica graphics

Or with Piecewise

Mean@TransformedDistribution[
  Piecewise[{{w + x1, w > c}, {w + x2, w < c}}],
  {
   Distributed[w, GammaDistribution[k0, th0]]
   , Distributed[x1, GammaDistribution[k1, th1]]
   , Distributed[x2, GammaDistribution[k2, th2]]
   }]
Piecewise[{{k0*th0 + k1*th1, c <= 0}}, ((k0*th0 + k2*th2)*Gamma[k0] + (k1*th1 - k2*th2)*Gamma[k0, c/th0])/Gamma[k0]]

Mathematica graphics

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}}],
    {
     Distributed[w, GammaDistribution[k0, th0]]
     , Distributed[x1, GammaDistribution[k1, th1]]
     , Distributed[x2, GammaDistribution[k2, th2]]
     }]
  ]]
(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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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