Utilities


Isserlis theorem from cumulant formulas

Step 1: Get expression for 4th moment in terms of cumulants
In[]:=
cumulants=convert[Moment[{1,1,1,1}]];​​cumulants//TraditionalForm
Out[]//TraditionalForm=
κ
0,0,0,1
κ
0,0,1,0
κ
0,1,0,0
κ
1,0,0,0
+
κ
0,0,1,1
κ
0,1,0,0
κ
1,0,0,0
+
κ
0,0,1,0
κ
0,1,0,1
κ
1,0,0,0
+
κ
0,0,0,1
κ
0,1,1,0
κ
1,0,0,0
+
κ
0,1,1,1
κ
1,0,0,0
+
κ
0,0,1,0
κ
0,1,0,0
κ
1,0,0,1
+
κ
0,1,1,0
κ
1,0,0,1
+
κ
0,0,0,1
κ
0,1,0,0
κ
1,0,1,0
+
κ
0,1,0,1
κ
1,0,1,0
+
κ
0,1,0,0
κ
1,0,1,1
+
κ
0,0,0,1
κ
0,0,1,0
κ
1,1,0,0
+
κ
0,0,1,1
κ
1,1,0,0
+
κ
0,0,1,0
κ
1,1,0,1
+
κ
0,0,0,1
κ
1,1,1,0
+
κ
1,1,1,1
Step 2: Set to zero all high-order cumulants
In[]:=
truncated=dropTerms[cumulants,rank[#]>2&];​​truncated//TraditionalForm
Out[]//TraditionalForm=
κ
0,0,0,1
κ
0,0,1,0
κ
0,1,0,0
κ
1,0,0,0
+
κ
0,0,1,1
κ
0,1,0,0
κ
1,0,0,0
+
κ
0,0,1,0
κ
0,1,0,1
κ
1,0,0,0
+
κ
0,0,0,1
κ
0,1,1,0
κ
1,0,0,0
+
κ
0,0,1,0
κ
0,1,0,0
κ
1,0,0,1
+
κ
0,1,1,0
κ
1,0,0,1
+
κ
0,0,0,1
κ
0,1,0,0
κ
1,0,1,0
+
κ
0,1,0,1
κ
1,0,1,0
+
κ
0,0,0,1
κ
0,0,1,0
κ
1,1,0,0
+
κ
0,0,1,1
κ
1,1,0,0
Step 3: Convert back to moments and simplify
In[]:=
moments=convertAll[truncated]//Simplify;​​moments//TraditionalForm
Out[]//TraditionalForm=
-2
μ
0,0,0,1
μ
0,0,1,0
μ
0,1,0,0
μ
1,0,0,0
+
μ
0,1,1,0
μ
1,0,0,1
+
μ
0,1,0,1
μ
1,0,1,0
+
μ
0,0,1,1
μ
1,1,0,0
Visualize
In[]:=
visualize[moments]
Out[]=
-2
+
+
+

Example for third moment

In[]:=
cumulants=convert[Moment[{1,1,1}]];​​truncated=dropTerms[cumulants,rank[#]>2&];​​moments=convertAll[truncated]//FullSimplify//Expand;​​moments//TraditionalForm​​visualize[moments]
Out[]//TraditionalForm=
-2
μ
0,0,1
μ
0,1,0
μ
1,0,0
+
μ
0,1,1
μ
1,0,0
+
μ
0,1,0
μ
1,0,1
+
μ
0,0,1
μ
1,1,0
Out[]=
-2
+
+
+

Convert 4 Gaussian einsums into matrix notation

In[]:=
Clear[A,X,d,mu]​​Assuming[(X|A)∈Matrices[{d,d}]&&mu∈Vectors[d],​​Print[FromTensor@TensorReduce@einsum[{{2,3},{1,2},{3,4}}{1,4},A,X,X]//TraditionalForm];​​Print[FromTensor@TensorReduce@einsum[{{2,3},{1,3},{2,4}}{1,4},A,X,X]//TraditionalForm];​​Print[FromTensor@TensorReduce@einsum[{{2,3},{1,4},{2,3}}{1,4},A,X,X]//TraditionalForm];​​Print[FromTensor@TensorReduce@einsum[{{2,3},{1},{2},{3},{4}}{1,4},A,mu,mu,mu,mu]//TraditionalForm];​​]​​
X.A.X
X.
T
A
.X
XTr[
T
A
.X]
mu.A.mumumu

Verify final formula

In[]:=
CircleTimes=KroneckerProduct;​​d=3;​​Σ=RandomInteger[{-5,5},{d,d}];​​Σ=Σ.Transpose[Σ]+IdentityMatrix[d];(*ensurepositivedefiniteness*)​​μ=RandomInteger[{-5,5},{d}];​​Assert[Min[N@Eigenvalues[Σ]]>0];​​dist=MultinormalDistribution[μ,Σ];​​Ex[expr_]:=Expectation[expr,xdist];(*Expectationwithrespecttox*)​​​​​​Clear[xx,aa];​​x=Array[xx,d];​​A=Array[aa,{d,d}];​​​​X2=Ex[x⊗x];(*uncenteredsecondmoment,sameasΣ+μ⊗μ*)​​​​Ex[((x⊗x).A.(x⊗x))]==X2.A.X2+X2.A.X2+X2Tr[A.X2]-2μ.A.μμ⊗μ//Reduce
Out[]=
True
In[]:=
deploy
Sun 29 Aug 2021 11:27:33