Cumulants in terms of Moments

This notebook implements my simple recurrence relation to generate expressions for cumulants in terms of power moments.
It was derived in  “Application of the Ratio Method of EXAFS Analysis to Disordered Systems”, 
Grant Bunker, Nuclear Instruments and Methods in Physics Research, (1983)
If you use it, please cite it.

writefile=False;(*flag whether to write out the expressions to a file in the same directory (folder) as this notebook*)

nmax=50;(*up to 50 takes a bit less than a minute on MacBook Pro M2 Max - see timing benchmarks later in this notebook*)

n is the order of the various moments 
γ is a formal parameter of no significance in itself  - it goes away
this employs an unnormalized distribution (upper case P[n,γ]);
expressions are generated, then related to moments ((Subscript[p, n])) of normalized distribution.

this calculates the expressions for all cumulants up to nmax in terms of power moments
the expressions grow rapidly with nmax
results are stored in “tableau”
the tableau can be extended to higher order simply by applying the nextn function to the last element of tableau and continuing it.
Each successive evaluation takes about 1.2274 times as long as the previous one (until memory limits are hit presumably)
This is the whole thing right here:

nmax=50;(*maximum order of cumulant expression to generate*)
(*each increase in nmax by 5 required about a factor of 3 in execution time*)
nextn=D[#,γ]/.{D[P[n_,γ],γ]->P[n+1,γ]}&;
Print["timing=",AbsoluteTiming[tableau=NestList[nextn,Log[P[0,γ]],nmax];]]

here is the symbolic representation that is suitable for continuing to higher order simply by applying “nextn”

tableau[[5]]

it works!

nextn@tableau[[5]]==tableau[[6]]

this expresses the tableau in more readable form (nicer forms are below)

cumulants=Drop[tableau//.{P[n_,γ]->p[n],p[0]->1},1];

cumulants[[6]]

nicer  yet

cumulants[[6]]/.p[n_]->(Subscript[p, n])

this  writes  out  the  file  to the directory in which this notebook is located - it  can  take  some  space, so be warned

If[writefile,SetDirectory[NotebookDirectory[ ]];Export["cumulants.txt",cumulants,"Lines"]

this verifies that the data are readable (and how to do so) and shows whether the re-read file agrees with the original

If[writefile,test=ToExpression@Import["cumulants.txt","Lines"]; cumulants==test]

the rest of this is purely informational

the number of terms grows rapidly with n;   c[50] has more than 200,000  terms

nterms=Length/@List@@@tableau;(*number of terms in each expression, vs n*)
(*grows a little slower than exponentially*)
GraphicsRow[{ListLogPlot[nterms,Frame->True,GridLines->Full],ListPlot[Ratios[nterms]]}]

Show[{ListLogPlot[nterms,Frame->True,GridLines->Full],LogPlot[Exp[fitfun/.soln],{n,1,50}]}]

each successive evaluation takes about 1.24 times as long as the one before

Through[{Mean,StandardDeviation}[Ratios@Table[AbsoluteTiming[nextn@tableau[[n]];][[1]],{n,30,40}]]]

ByteCount[tableau]/1024.^3 GB

MaxMemoryUsed[]/1024.^3 GB

here are some nicely formatted results

p[n] is the nth moment (arbitrary origin)
here we limit the order so as  to fit on the page

cumulants[[1;;6]]/.{p[n_]->(Subscript[p, n])}//ColumnForm

this does it in terms of central moments (those calculated with origin at the centroid)
here we limit the printout  to fit on the page nicely
but the calculated ones here go up to n=50 and easily can go somewhat higher subject to time and memory limitations

format the result little more nicely

nicetable=Transpose[{Range[1,7],cumulants[[1;;7]]/.{p[n_]->(Subscript[p, n])}}];

the  following  expresses  cumulants  in  terms  of  moments  relative  to  arbitrary  origin

nicetable//TableForm

here are expressions for cumulants in terms of central moments

nicetable/.(Subscript[p, 1])->0//TableForm

some more statistics on complexity and execution times

the complexity of the expressions grows rapidly with m aka nmax
 ln(leafcount) itself grows as a power law with the order of the cumulant 
 the recursion relation generates the expressions for all lower orders if you want them

fig=ListLogLogPlot[{Log[LeafCount/@tableau]}]

soln=FindFit[Log[LeafCount/@tableau],a lnm^b,{a,b},lnm]

AbsoluteTiming[Nest[nextn,Log[P[0,γ]],35]//.{P[n_,γ]->p[n],p[0]->1};]

timetocomputeupto[m_]:=AbsoluteTiming[NestList[nextn,Log[P[0,γ]],m]//.{P[n_,γ]->p[n],p[0]->1};][[1]]

timetocomputeupto[35]

timetocomputeupto[40]

%/%%(*increasing n by 5 increases the compute time by about a factor of 3, 1.24^5*)

MemoryInUse[]/1024.^3 GB

MaxMemoryUsed[]/1024.^3 GB

ByteCount[tableau]/1024.^3 GB

ByteCount[cumulants]/1024.^3 GB

LeafCount[tableau]

LeafCount[cumulants]