​
In[]:=
Clear["Global`*"];​​Clear[b,A,n,p,α,Λ,X,m1,pp,mm,a,XL,λs,tra,frob,varyingp,x];​​​​n=2;(*here,nshouldbegreaterthanorequaltotheactualnumberofeigenvalueswewant*)​​​​XL[c_,λ_,n_]:=Table[PDF[c,x/λ[[i]]]/Abs[λ[[i]]],{i,1,n}]//FullSimplify(*XLmakesalistofdistributionsoftherandomvariablesx_iλ_i*)​​Λ=Refine[Table[Subscript[λ,i],{i,n}],For[i=1,i<n+1,i++,Subscript[λ,i]∈]];​​A=DiagonalMatrix[Λ];(*WecanrunnormfunctionsforanyHermitianmatriceswewant,buthereΛisjustalistofλ_isandAitsdiagonalmatrix*)​​​​$Assumptions=Λ[[1]]≠0&&Λ[[2]]≠0&&Λ[[3]]≠0&&Λ[[4]]≠0&&e<f&&q<s&&s<r;​​​​(*zetamgf[X_,A_,n_]=Product[Zeta[λs[A][[i]]t],{i,1,n}];*)​​mgfa[X_,A_,n_]:=Product[MomentGeneratingFunction[If[ListQ[X],X[[i]],X],Eigenvalues[A][[i]]t],{i,1,n}];(*Producesthemomentgeneratingfunctionasdescribedinthepaper.*)​​​​Series[mgfa[X,A,n],{t,0,3}];(*ProvidesaTaylorserieswith3terms.Ifyouwantfewer(probably2)eigenvalues,youcanchangetheinputfromnto2,orredefinenabove.*)​​CoefficientList[%,t]//FullSimplify​​normp[X_,A_,n_,p_]:=(Expectation[Abs[Sum[Subscript[B,i],{i,1,n}]]^p,Table[Subscript[B,i]ProbabilityDistribution[XL[If[ListQ[X],X[[i]],X],Eigenvalues[A],n][[i]],{x,-Infinity,Infinity}],{i,1,n}]]/p!)(*GivesthepthnormofA*)​​​​normp[NormalDistribution[2,1],A,2,2]​​normp[NormalDistribution[2,1],A,2,4]​​normp[NormalDistribution[2,1],A,2,6]​​​​normp[NormalDistribution[0,1],A,2,2]​​normp[NormalDistribution[0,1],A,2,4]​​normp[NormalDistribution[0,1],A,2,6]​​​​​​hardxyed[X_,A_,n_,p_]:=normp[X,A,n,p]/.
λ
1
z/.
λ
2
x/.
λ
3
y​​xyed[X_,A_,n_,p_]:=SeriesCoefficient[Series[mgfa[X,A,n],{t,0,p}],p]/.
λ
1
z/.
λ
2
x/.
λ
3
y​​​​xyer[X_,A_,n_,p_]:=If[EvenQ[p],xyed[X,A,n,p],hardxyed[X,A,n,p]]​​xyedA2[X_,p_]:=xyer[X,A,2,p]​​xyedA3[{X_,p_}]:=xyer[X,A,3,p]​​absoluteunit[f_]:=Solve[f1,z]​​​​parameterize[f_,p_]:=(f/.xCos[t]/.zSin[t])^(1/p)​​parameterize3d[f_,p_]:=(f/.xCos[t]Sin[u]/.zSin[t]Sin[u]/.yCos[u])^(1/p)​​axes2d[x_]:={Cos[t]/x,Sin[t]/x}​​axes3d[x_]:={Cos[t]Sin[u]/x,Sin[t]Sin[u]/x,Cos[u]/x}​​​​paraplot[f_,p_]:=ParametricPlot[axes2d[parameterize[f,p]],{t,0,2Pi}]​​multiparaplot[ff_,pp_,ll_]:=ParametricPlot[Evaluate[Transpose[axes2d[MapThread[parameterize,{ff,pp}]]]],{t,0,2Pi},PlotLabelsll,AxesLabel{
x
1
,
x
2
}]​​paraplot3d[f_,p_]:=ParametricPlot3D[axes3d[parameterize3d[f,p]],{t,0,2Pi},{u,-Pi,Pi}]​​ultimateplot2d[X_,p_]:=paraplot[xyedA2[X,p],p]​​multiultimateplot[XX_,pp_,ll_]:=multiparaplot[MapThread[xyedA2,{XX,pp}],pp,ll]​​ultimateplot3d[X_,p_]:=paraplot3d[xyedA3[X,p],p]​​label[lab_,list_]:=Table[StringJoin[{lab," = ",ToString[list[[i]]]," "}],{i,Length[list]}]​​​​(*zetaxyed[X_,A_,n_,p_]:=SeriesCoefficient[Series[zetamgf[X,A,n],{t,0,p}],p]/.
λ
1
z/.
λ
2
x/.
λ
3
y​​zetaxyedA2[X_,p_]:=zetaxyer[X,A,2,p]​​zetamultiultimateplot[XX_,pp_,ll_]:=multiparaplot[MapThread[zetaxyedA2,{XX,pp}],pp,ll]​​Series[zetamgf[X,A,n],{t,0,3}]*)​​​​
Part
:Part 3 of {
λ
1
,
λ
2
} does not exist.
Part
:Part 4 of {
λ
1
,
λ
2
} does not exist.
Out[]=

2
MomentGeneratingFunction[X,0]
,MomentGeneratingFunction[X,0](
λ
1
+
λ
2
)
(0,1)
MomentGeneratingFunction
[X,0],
λ
1
λ
2
2
(0,1)
MomentGeneratingFunction
[X,0]
+
1
2
MomentGeneratingFunction[X,0](
2
λ
1
+
2
λ
2
)
(0,2)
MomentGeneratingFunction
[X,0],
1
6
3
λ
1
λ
2
(
λ
1
+
λ
2
)
(0,1)
MomentGeneratingFunction
[X,0]
(0,2)
MomentGeneratingFunction
[X,0]+MomentGeneratingFunction[X,0](
3
λ
1
+
3
λ
2
)
(0,3)
MomentGeneratingFunction
[X,0]
Out[]=
5
2
λ
1
+8
λ
1
λ
2
+5
2
λ
2
2Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Out[]=
43
4
λ
1
+112
3
λ
1
λ
2
+150
2
λ
1
2
λ
2
+112
λ
1
3
λ
2
+43
4
λ
2
24Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Out[]=
499
6
λ
1
+1704
5
λ
1
λ
2
+3225
4
λ
1
2
λ
2
+3920
3
λ
1
3
λ
2
+3225
2
λ
1
4
λ
2
+1704
λ
1
5
λ
2
+499
6
λ
2
720Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Out[]=
2
λ
1
+
2
λ
2
2Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Out[]=
2
(
2
λ
1
+
2
λ
2
)
8Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Out[]=
3
(
2
λ
1
+
2
λ
2
)
48Abs[
λ
1
]Abs[
λ
2
]
1
2
λ
1
1
2
λ
2
Syntax
:"(*zetaxyed[X_,A_,n_,p_]:=SeriesCoefficient[Series[zetamgf[X,A,n],{t,0,p}],p]/.
λ
1
z/.
λ
2
x/.
λ
3
yzetaxyedA2[X_,p_]:=zetaxyer[X,A,2,p]1Series[zetamgf[X,A,n],{t,0,3}]*) is incomplete; more input is needed.
alphasToTry={2.1,3,4,10};​​musToTry={-2,-1,0,1,6};​​psToTry={2,4};​​bernoulliParametersToTry={.01,.5,.99};​​​​fixedP=10;​​fixedBernoulliParameter=.5;​​fixedAlpha=4.1;​​​​paretosVaryingAlpha={Table[ParetoDistribution[1,alphasToTry[[i]]],{i,Length[alphasToTry]}],Table[fixedP,Length[alphasToTry]],CreateLabels["α",alphasToTry]};​​paretosVaryingP={Table[ParetoDistribution[1,fixedAlpha],Length[psToTry]],psToTry,CreateLabels["p",psToTry]}​​exponentialsVaryingP={Table[ExponentialDistribution[1],Length[psToTry]],psToTry,CreateLabels["p",psToTry]}​​normalsVaryingMean={Table[NormalDistribution[musToTry[[i]],1],{i,Length[musToTry]}],Table[fixedP,Length[musToTry]],label["μ",musToTry]};​​normalsVaryingP={Table[NormalDistribution[0,1],Length[psToTry]],psToTry,CreateLabels["p",psToTry]};​​BernoullisVaryingParameter={Table[BernoulliDistribution[bernoulliParametersToTry[[i]]],{i,Length[bernoulliParametersToTry]}],Table[fixedP,Length[bernoulliParametersToTry]],bernoulliParametersToTry};​​BernoullisVaryingP={Table[BernoulliDistribution[fixedBernoulliParameter],Length[psToTry]],psToTry,psToTry};​​schattenPNorms={Table[Abs[x]^psToTry[[i]]+Abs[z]^psToTry[[i]],{i,Length[psToTry]}],psToTry,CreateLabels["x",{1}]}​​schattenInfinityNorm={{Max[Abs[x],Abs[z]]},{1},{"p = ∞"}}​​​​blob={{{NormalDistribution[-2,1],BernoulliDistribution[.1]}},{10},{StringJoin[ToString[Subscript[x,1]+Subscript[x,2],StandardForm]," = 10"]}}​​scaledpnorms={Table[Abs[x/7.1775512683]^otherpsToTry[[i]]+Abs[z/7.1775512683]^otherpsToTry[[i]],{i,Length[otherpsToTry]}],otherpsToTry,label["p",otherpsToTry]}​​mix[Y_,X_]:=Join[Y,{MapThread[xyedA2,{X[[1]],X[[2]]}],X[[2]],X[[3]]},2]​​pmix=Join[pnorms,infinitynorm,2];​​blobmix=mix[straightline,blob];​​​​​​multiultimateplot@@BernoullisVaryingP​​​​
Out[]=
{{ParetoDistribution[1,4.1],ParetoDistribution[1,4.1]},{2,4},{p = 2 ,p = 4 }}
◼
  • ​