Symbolic Mean and Variance Calculations

Question

I need a tutorial that gets me started on basic symbolic calculations with random variables. My naive attempts did not get me far.

In[1]:= x= NormalDistribution[m,s];

In[2]:= y= NormalDistribution[m,s];

In[3]:= Mean[x+y] // FullSimplify
Out[3]:= Mean[2NormalDistribution[m,s]]

I would have expected to see 2*m here as a result, but this was obviously my fail.

My goal is to be able to evaluate expressions such as Var[ w_i * R_i ] and get back Sum_i[ Sum_j[ 2 * w_i * w_j * Cov[R_i,R_j] ]]. This means I have to tell Mma that w_i are constant and R_i are random variables. (In this case, I wouldn't have to tell that the variables are normally distributed, either; just that they are RVs.) Ideally, I would subsequently be able to tell Mma that Cov[R_i,R_j]=0 for i != j and get back Sum_i w_i Var[ R_i ].

because I often work with repeated draws, being able to specify that all R_i are iid draws from the same distribution would be helpful, too.

Pointers would be highly appreciated. (As would an example...)

Answer 1

If you're only interested in the means and variances of linear combinations of random variables (where the means, variances, and covariances exist), then consider the following.

Suppose the means and covariances are $\mu$ $\Sigma$, respectively, then the linear combination of the multivariate random variable $X$ and another known vector $\omega$ of the same dimension has mean $$\omega.\mu$$ and the variance is $\omega'.\Sigma.\omega$. Using Mathematica notation:

n = 3;
μ = Table[m[i], {i, n}];
Σ = 
  Table[cov[Min[i, j], Max[i, j]], {i, n}, {j, n}] /. 
   cov[i_, i_] -> σ[i]^2;
ω = Table[w[i], {i, n}];

(* Mean of ω.X *)
mean = ω . μ
(* m[1] w[1]+m[2] w[2]+m[3] w[3] *)

(* Variance of ω.X *)
var = ω . Σ . ω // Simplify 
(* 2 cov[1,2] w[1] w[2]+2 cov[1,3] w[1] w[3]+2 cov[2,3] w[2] \
w[3]+w[1]^2 σ[1]^2+w[2]^2 σ[2]^2+w[3]^2 σ[3]^2 *)

If you're interested in more complicated functions, please include an example in your question.

Addition:

Your example in a comment:

$$Z=\prod _{i=1}^3 (a+b X_i+\epsilon_i)$$

To use replacement rules to determine $E(Z)$ and $Var(Z)$ one could construct the following:

z = Product[a + b x[i] + e[i], {i, 3}] // Expand
(* a^3 + a^2 e[1] + a^2 e[2] + a e[1] e[2] + a^2 e[3] + a e[1] e[3] + 
 a e[2] e[3] + e[1] e[2] e[3] + a^2 b x[1] + a b e[2] x[1] + 
 a b e[3] x[1] + b e[2] e[3] x[1] + a^2 b x[2] + a b e[1] x[2] + 
 a b e[3] x[2] + b e[1] e[3] x[2] + a b^2 x[1] x[2] + 
 b^2 e[3] x[1] x[2] + a^2 b x[3] + a b e[1] x[3] + a b e[2] x[3] + 
 b e[1] e[2] x[3] + a b^2 x[1] x[3] + b^2 e[2] x[1] x[3] + 
 a b^2 x[2] x[3] + b^2 e[1] x[2] x[3] + b^3 x[1] x[2] x[3] *)

Now assuming that the random variables are all independent of each other the mean of $Z$ is

rules = {e[i_]^2 -> σ^2, e[i_] -> 0, x[i_]^2 -> σ[i]^2 + μ[i]^2, x[i_] -> μ[i]};
ez = z /.rules // FullSimplify
(* (a + b μ[1]) (a + b μ[2]) (a + b μ[3]) *)

The variance is

vz = ((z^2 // Expand) /. rules) - ez^2 // FullSimplify
(* σ^6 + b^2 σ^4 (μ[1]^2 + μ[2]^2 + μ[3]^2 + σ[1]^2 + σ[2]^2 + σ[3]^2) +
   b^4 σ^2 (μ[3]^2 σ[1]^2 + μ[3]^2 σ[2]^2 + σ[1]^2 σ[2]^2 + (σ[1]^2 + σ[2]^2) σ[3]^2 +
   μ[2]^2 (μ[3]^2 + σ[1]^2 + σ[3]^2) + μ[1]^2 (μ[2]^2 + μ[3]^2 + σ[2]^2 + σ[ 3]^2)) + 
   a^4 (3 σ^2 + b^2 (σ[1]^2 + σ[2]^2 + σ[3]^2)) + b^6 ((μ[1]^2 + σ[1]^2) σ[2]^2 (μ[3]^2 + σ[3]^2) +
   μ[2]^2 (μ[3]^2 σ[1]^2 + (μ[1]^2 + σ[1]^2) σ[3]^2)) + 
   2 a^3 (2 b σ^2 (μ[1] + μ[2] + μ[3]) + 
   b^3 (μ[3] (σ[1]^2 + σ[2]^2) + μ[2] (σ[1]^2 + σ[3]^2) + μ[1] (σ[2]^2 + σ[3]^2))) + 
   2 a (b σ^4 (μ[1] + μ[2] + μ[3]) + b^5 (μ[3] (μ[2] (μ[2] + μ[3]) σ[1]^2 + (μ[1] (μ[1] + μ[3]) + 
   σ[1]^2) σ[2]^2) + (μ[2] (μ[1] (μ[1] + μ[2]) + σ[1]^2) + 
   μ[1] σ[2]^2) σ[3]^2) + b^3 σ^2 (μ[2]^2 μ[3] + μ[1]^2 (μ[2] + μ[3]) +
   μ[3] (σ[1]^2 + σ[2]^2) + μ[2] (μ[3]^2 + σ[1]^2 + σ[3]^2) + 
   μ[1] (μ[2]^2 + μ[3]^2 + σ[2]^2 + σ[3]^2))) + 
   a^2 (3 σ^4 + 2 b^2 σ^2 ((μ[1] + μ[2] + μ[3])^2 + σ[1]^2 + σ[2]^2 + σ[3]^2) + 
   b^4 (4 μ[1] μ[3] σ[2]^2 + (μ[1]^2 + σ[1]^2) σ[2]^2 + μ[3]^2 (σ[1]^2 + σ[2]^2) + 
   (μ[1]^2 + σ[1]^2 + σ[2]^2) σ[3]^2 + μ[2]^2 (σ[1]^2 + σ[3]^2) + 
   4 μ[2] (μ[3] σ[1]^2 + μ[1] σ[3]^2))) *)

The rules would need to account for any covariances or higher order joint moments which you say you'll add after the expansion.