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]]
  • $\begingroup$ Your function "invert" does not use its argument. $\endgroup$ Mar 20, 2023 at 9:35
  • $\begingroup$ Also I think your m[t] is not correct for defining the tth moment. I think you want the sum of i^t*h[i] for i in the range of the distribution. $\endgroup$ Mar 20, 2023 at 16:39
  • $\begingroup$ I may be using the wrong terminology. $m(k)$ gives $E_X[X^0 p(X)^k]$ rather than $E_X[X^k]$ (is there a name for it?). It seems there's a formula to invert the mapping from $m(1),\ldots,m(d)$ although it's not likely to work for the example above $\endgroup$ Mar 20, 2023 at 17:03

1 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]= \
308091434897742274739090092408039663230100469580*10^-117 *)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.