Symbolic Mean and Variance Calculations

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...)

• look up TransformedDistribution in the docs. E.g., try Mean[TransformedDistribution[a + b, {Distributed[a, x], Distributed[b, y]}]]
– kglr
Commented May 30, 2021 at 20:06
• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Commented May 30, 2021 at 21:16
• great stuff in the TransformedDistribution link (thank you), but not really getting me to where I need to get to---which is manipulation of means and variances. :-( Commented May 31, 2021 at 3:04
• You need to keep the notations for a random variable and a distribution separate. As such 2*NormalDistribution[m,s] is not a legitimate command. If you want to avoid TransformedDistribution, consider using Expectation.
– JimB
Commented May 31, 2021 at 5:01

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.

$$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.

• Mille grazie. I need Mma to work with RVs. indeed, my example was too simple. let's say I wanted to calculate Var[ \Prod_i[ (a + b*x[i] + e[i]) ] ], and a and b are constants while x[i] and e[i] are RVs . there will be tons of cross-terms, and Mma will help me not lose them. I then will want to add assumptions, like all e[i] are iid errors. Commented May 31, 2021 at 16:30
• Maybe using replacement rules after expanding to all of the cross-product terms. I'll add in your product example but I think you're going to have to get used to how Mathematica does things (just like any other language). I would love to have it deal directly with random variables (without having to specify a specific distribution). But there are folks in this forum far better than me. So maybe there's hope.
– JimB
Commented May 31, 2021 at 18:53
• Thank you. This goes a very long way towards where I want to go. It's interesting that Mathematica does not have it "built-in". Commented Jun 2, 2021 at 14:56
• @ivoWelch Do you know of any other software that does have it "built-in" ? I'd like to use that software.
– JimB
Commented Jun 2, 2021 at 14:58
• me, too. no, I do not know of any other ones. incidentally, you have some syntax errors in the rules example. Commented Jun 2, 2021 at 18:02