Moment problem for a discrete probability distribution

Question

Is there an elegant way to use Mathematica to find entries of a discrete probability distribution $\{p_i\}$ by using a finite number of calls to $m(t)=\sum_i p_i^t$?

There's a wikipedia article on the moment problem but it's fairly abstract.

IE, for example, for the code below, successful invert[m] would output h

h = 1./Range[10000];
h = h/Total[h];
m[t_] := Total[h^t];

invert[m_] := (
   ConstantArray[1, Length[h]]
   );

Print["Error is ", Norm[invert[m] - h]]

Answer 1

In principle this can be inverted by solving a linear system. In practice it might be difficult for the case of a large domain, and I would expect one might need a black-box iterative solver.

Here is the same example except with probability measure defined on Range[5]. Note that I am using a corrected moment computation below.

max = 5;
h = 1./Range[max];
h = h/Total[h];
m[n_] := h . (Range[max]^n);

First let's check the actual probability distribution.

In[740]:= Table[h[[j]], {j, max}]

(* Out[740]= {0.437956, 0.218978, 0.145985, 0.109489, 0.0875912} *)

The right-hand-side of the linear system will be the moments. I will use the zero through fourth because it makes the matrix construction slightly simpler.

rhs = Table[m[j], {j, 0, max - 1}]

(* Out[749]= {1., 2.18978, 6.56934, 24.0876, 98.5401} *)

A bit of thought shows that the matrix we need in order to set up the appropriate system is the transpose of a Vendermonde matrix defined on the support of the probability measure. We will next use that to solve the system.

vmat = VandermondeMatrix[Range[max]];
pj = LinearSolve[Transpose@vmat, rhs]

(* Out[751]= {0.437956, 0.218978, 0.145985, 0.109489, 0.0875912} *)

Note that we have in fact recovered the probability distribution.

Some testing indicates that we need to bump the precision to scale with dimension if we want reliable results. This of course means the solve time will rise substantially with dimension. Here is what I get at max=1000. I assess error by the maximum component of the residual.

m[n_, max_] := h . (Range[max]^n);
max = 1000;
h = 1/Range[max];
h = N[h/Total[h], max];

rhs = Table[m[j, max], {j, 0, max - 1}];
vmat = VandermondeMatrix[Range[max]];
Timing[pj = LinearSolve[Transpose@vmat, rhs];]
Max[Abs[pj - Table[h[[j]], {j, max}]]]

(* Out[813]= {6.96261, Null}

Out[814]= \
5.54682699290347132366334935835377862569126864261898930592167088339721\
6535874428170235873486394904117700675166858326547221444890544975590619\
5592941877951826436679154857249440225552646668109589531131194083595479\
6180460454984157584498672271023556602999130161221006457906685102016423\
1746743679648741153409418099382190627469397083934113609253526505366745\
9812595544941701077190477304212378059202466992076424069385724859593883\
0964208138994311146419785285562138374443922728104683934917887149620746\
6140575547300389533048381282608195517041685477069237325500653906104987\
9561366442202883775766130040759073477328163091560136111691769172306175\
8469259487156068204542442459360764924597159847753622838936239374468685\
2565974882699291981983509947770267917952400799677091557733748690067802\
2209540697710929780857463084268457402190894496997545921534643942945935\
308091434897742274739090092408039663230100469580*10^-117 *)