In[]:=
Expectation[m^2*n*o,{m,n,o}MultinormalDistribution[{u1,u2,u3},{{s1^2,s12,s13},{s12,s2^2,s23},{s13,s23,s3^2}}]]
Out[]=
2
μ
1
μ
2
μ
3
+
μ
2
μ
3
2
σ
1
+2
μ
1
μ
3
σ
12
+2
μ
1
μ
2
σ
13
+2
σ
12
σ
13
+
2
μ
1
σ
23
+
2
σ
1
σ
23
​