How to estimate the probability mass function of a discrete variable from moments

Question

Consider a bounded, discrete random variable $X$ whose range is $(0,1,\ldots, M)$. We are given the first $k$ moments of this distribution, call them $m_1, \ldots, m_k$. We are interested in estimating the probability mass function of X using this information.

A simple, first principle approach for this seems to be to write down a linear equation for each $m_j$, and solve the system of linear equations. Let $A_{k \times M}$ denote the matrix whose $j^{th}$ row is $(1, 2^j, \ldots, M^j)$, let $p$ be the PMF vector of length $M$, i.e., $p_i=Prob[X=i]$, and $m$ be the vector $(m_1, \ldots, m_k)$. Then we have to solve the linear system $Ap = m$. Of course, in general this won't have an unique solution, but let's say we are happy with any valid $p$ that solves this system.

My question is: how can one enforce the constraints $p_i \geq 0$ for all $i$, and $\sum p_i = 1$? For example, is there an easy way to use linear programming type ideas here?

As a follow up, what are the pros and cons of this approach compared to using the moments to construct an (imperfect) moment generating function, and finding the PMF from the MGF?

Answer 1

You are seeking to solve the underdetermined system of linear equations $\mathbf{A} \mathbf{p} = \mathbf{m}$ subject to the additional constraint that $\mathbf{p}$ must be a probability vector (i.e., $p_i \geqslant 0$ and $\sum p_i = 1$). Finding a useful characterisation of the full set of solutions is quite a difficult programming problem, but it should be fairly simple to find a single solution in a particular case. For the latter problem, assuming there is a solution that is not at a boundary point (which will usually be the case), all you really need to do is to get to a point that is "close" to a solution point that obeys the required constraints and then use this to fix $M-k+1$ of the values in the system, and then solve the remaining values using a standard fully-determined system of linear equations.

There are many ways to solve this problem. One way is to rewrite it as an unconstrained non-linear problem, and iterate towards a solution. To do this, define the unconstrained vector $\mathbf{r} = (r_1,...,r_M) \in \mathbb{R}^M$ (and let $r_0 \equiv 0$) and use the softmax transformation to obtain the associated probability vector:

$$p_i = \frac{\exp(r_i)}{1+\sum_{j=0}^M \exp(r_j)} \quad \quad \quad \text{for } i = 0,...,M.$$

You can then write your moment equations as:

$$m_a = \frac{\sum_{i=1}^M i^a \cdot \exp(r_i)}{1+\sum_{i=1}^M \exp(r_i)} \quad \quad \quad \quad \quad \text{for } a = 1...,k. $$

You can use non-linear programming to iterate towards a solution to this system of non-linear equations, and thereby obtain a point $\hat{\mathbf{r}}$ that is close to a solution point. Create the corresponding probability vector $\hat{\mathbf{p}}$ and take the values $\hat{p}_k,...,\hat{p}_M$ from this vector. You now have a fully-determined system of linear equations $\mathbf{A} \mathbf{p} = \mathbf{m}$ where the values $\hat{p}_k,...,\hat{p}_M$ are fixed.

Assuming that there is a solution to your system that is not on the boundary point, it should be the case that there will be a solution that is "close enough" to $\hat{\mathbf{p}}$ so that holding its last values constant gives a solution point from the corresponding system of linear equations. If you add some specific numbers to your problem then it should be fairly simple to find a solution by this technique.

Answer 2

I have worked with this problem a lot, and I cannot find a satisfactory answer, but here are a few ways that I've attacked it.

(1) If $M=k$, and if you know that a valid probability distribution can be created, then this will be the solution: $$ \begin{equation} \begin{bmatrix} 1 & ... & 1\\ 0 & ... & M\\ (0-mean)^2 & ... & (M-mean)^2\\ ...& ... & ... \end{bmatrix} P_{k} = \begin{bmatrix} 1 \\ mean \\ variance \\ .. \end{bmatrix} \end{equation} $$

(2) If $M>k$, then choose a subset of $[0,...,M]$ that is size $k$ that you know can be a domain for a valid probability distribution, and repeat (1).

For both (1) and (2), knowing you have domain that can produce a valid distribution isn't an easy task. Yamada and Primbs (in Construction of Multinomial Lattice Random Walks for Optimal Hedges) wrote some formulas to verify whether the domains are valid for moments up to kurtosis, but the formulas get pretty nasty for orders higher than that.

(3) As others have pointed out, this isn't a linear optimization problem: but you can turn it into a quadratic one. I use the Goldfarb-Idnani algorithm to minimize $pAp^T$, where $A$ is the identity, then add the constraints that $p_i\ge0$, $\sum_{i=0} ^{M} p_i = 1$, $\sum_{i=0} ^{M} ip_i = mean$, &etc. Goldfarb-Idnani doesn't require a guess that satisfies the constraints (which is your problem to begin with). It also minimizes entropy, which makes for a really nice looking solution. The downside is that it takes considerably longer than (1) or (2).

(4) Donghui Chen and Robert J.Plemmons, Nonnegativity Constraints in Numerical Analysis. Construct a matrix like in (1), but where $M>k$, then solve with their algorithm. Much faster than (3), but the distribution it gives will not be minimizing entropy. Slower than (1), but can work much better in practice.