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

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?

• Usually one solves this by minimizing some objective function, such as the negative entropy, subject to the constraints. But what are you asking? How to solve constrained linear problems? The question as you pose it is a linear program, so it's strange to see you ask whether "linear programming type ideas" would be applicable! – whuber Mar 27 '19 at 22:25
• My question is how to set up the constrained linear problem. For example, should one use the objective function $||Ap-m||+pen(\sum p_i - 1)$ where the second term penalizes for solutions that are not valid PMFs? How would this approach compare to using the negative entropy objective function?A follow-up question is, how does this optimization approach compare with finding the PMF from the moment generating function? – lordcretin Mar 28 '19 at 16:57
• Again, no penalty is needed: this is purely a linear program. Specifically, it's the first step of a linear program: find the set of feasible solutions. Your question doesn't describe any optimization problem at all. – whuber Mar 28 '19 at 18:04

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.