Wolfram Cloud Document

In[]:=

deploy

Thu 22 Sep 2022 18:23:10

Helper utilities

In[]:=

On[Assert];(* deploys with canonical name *)deploy:=Module[{notebookFn, parentDir,cloudFn,result},Print[DateString[]];notebookFn=FileNameSplit[NotebookFileName[]][[-1]];parentDir=FileNameSplit[NotebookFileName[]][[-2]];cloudFn=parentDir~StringJoin~"/"~StringJoin~notebookFn;result=CloudDeploy[SelectedNotebook[],CloudObject[cloudFn],Permissions"Public",SourceLinkNone];Print["Uploading to ",cloudFn];result];Unprotect[Ex];Clear[x,y,X,xx,yy,xx1,xx2];symmetricMatrix[d_,sym_]:=(Clear[sym];Array[sym[Min[#1,#2],Max[#1,#2]]&,{d,d}]);arbitaryMatrix[d_,sym_]:=(Clear[sym];Array[sym,{d,d}]);subVars[p_]:=With[{vars=Variables[p]},Thread[varsRandomInteger[{-5,5},Length[vars]]]]Ex[expr_]:=Expectation[expr,xdist]; (* Expectation with respect to x *)Exy[expr_]:=Expectation[expr,{xdist,ydist}]; (* Expectation with respect to x *)EX[expr_]:=Expectation[expr,{x1dist,x2dist}];Protect[Ex];computeDerived:=Module[{},Clear[aa,hh];Clear[xx];Assert[Head[dist]~endsWith~"Distribution"];d=Length@RandomVariate[dist];(* symbolic variables for vector and matrix X *);x=Array[xx,d];x1=Array[xx1,d];x2=Array[xx2,d];X={x1,x2};A=arbitaryMatrix[d,aa];As=symmetricMatrix[d,aa];H=arbitaryMatrix[d,hh];Hs=symmetricMatrix[d,hh];a=vec[A];as=vec[As];ii=IdentityMatrix[d];ones=ConstantArray[1,d];vars=Variables[A]~Join~Variables[H];X2=Ex[x⊗x];X2X2=Ex[(x⊗x).(x⊗x)];X4=Ex[Outer[Times,x,x,x,x]];X4flat=Flatten[X4,{{1,2},{3,4}}];Assert[X4flatEx[(x⊗x)⊗(x⊗x)]];];(* Setup symbolic and numeric matrices for checks *)subst:=(SeedRandom[1];Thread[varsRandomReal[{-1,1},{Length[vars]}]]);substInt:=(SeedRandom[1];Thread[varsRandomInteger[{-3,3},{Length[vars]}]]);(* setup distribution *)setupDistEmpirical:=(d=2;X=RandomInteger[{-5,5},{d,20}];dist=EmpiricalDistribution[Transpose@X];computeDerived;);setupDistGaussian:=(d=3;Σ=RandomInteger[{-5,5},{d,d}];Σ=Σ.Transpose[Σ]+IdentityMatrix[d];μ=RandomInteger[{-5,5},{d}];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[μ,Σ];computeDerived);setupDistCenteredGaussian:=(d=3;Σ=RandomInteger[{-5,5},{d,d}];Σ=Σ.Transpose[Σ]+IdentityMatrix[d];μ=ConstantArray[0,{d}];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[Σ];computeDerived);setupDistDiagonalGaussian:=(d=4;Σ=DiagonalMatrix[Table[1/i,{i,1,d}]];μ=Table[0,d];Assert[Min[N@Eigenvalues[Σ]]>0];dist=MultinormalDistribution[Σ];computeDerived)endsWith[a_,b_]:=StringMatchQ[ToString[a],___~~ToString[b]~~EndOfString];outer2[x_]:=Outer[Times,x,x];CircleTimes=KroneckerProduct;

Checking Gaussian formula

Centered Gaussian, m=1

In[]:=

SeedRandom[1];setupDistCenteredGaussian;mu2=outer2[μ];Clear[m];Print["symmetric A, centered Gaussian, m=1, stats.SE formula"];A=As;m=1;predicted=(m(m-1)Σ.A.Σ+m(m-1)(Σ.A.mu2+mu2.A.Σ)+m(m-1)μ.A.μmu2+mΣ.(A+A).Σ+m(Σ.(A+A).mu2+mu2.(A+A).Σ)+m(Tr[Σ.A]ii+μ.A.μ).(Σ+mu2));predicted=predicted/.substInt;Print[StringForm["Σ=``, μ=``, A=``,m=``",MatrixForm[Σ],MatrixForm[μ],MatrixForm[A/.substInt],m]];observed=Ex[outer2[x].A.outer2[x]]/.substInt;Print["Predicted: ",predicted//MatrixForm];Print["Observed: ",observed//MatrixForm];Print["matches: ",predicted==observed];

Helper utilities

Checking Gaussian formula

Final formula check

Symmetric formula check