# Decomposing Sums of Random Variables

Suppose I have $$M$$ random variables, and a number of realizations of each variable. Each RV has the probability mass function: $$\rho_{X_i}(x) = \begin{cases} p_i, & x = 1\\ 1-p_i, & x = 0\\ 0, & \text{otherwise} \end{cases}$$

Now consider $$Z$$, which is the sum of all $$X_i$$: $$Z = \sum_{i}X_i$$ How can I take $$Z$$ and find all $$p_i$$?

edit: A further note, I'd like to consider the case where you do not have information about how many $$X_i$$ there are.

### How I'm going about it:

Here is a picture of 30 samples of the time series $$Z_N(t)$$: $Z_n$(t)">

I have considered that, if $$\vec{Z_N}$$ is a vector which contains $$N$$ realizations of $$Z$$, then we can perhaps treat that as a time series, $$Z_N(t)$$, and perform an autocorrelation or a fourier transform, or something similar. I think that the expected wait time for $$X_i$$ ought to be $$\frac{1}{p_i}$$ samples, so I'd expected the ACF to have spikes at those intervals, but when I ran it, I got the following (this is with $$0 < p_i < 0.5$$):

• I'd probably take the generating function of the distribution of $Z$, and try to factor it into linear terms. Ideally, they should be of the form $1-p_i+p_iz$, but it might take a lot of samples to get to where that factorization becomes evident. – Brian Tung Mar 12 at 19:53
• I don't see why you would expect this time series to show any interesting autocorrelation. – Brian Tung Mar 12 at 19:56
• @BrianTung, Thanks for responding. Clearly the ACF is not meaningful, I was just showing how I was approaching it. I'm not very good at this kind of thing, but I figured that if the expected wait time between pulses from a given RV was 5, then you would expect to see a higher autocorrelation at 5 than at 4 or 6. For your suggestion, are you saying that we construct the PMF of $Z$ by taking many samples? – Fred Frey Mar 12 at 20:02
• The problem is that it's memoryless, so you won't see this. – Brian Tung Mar 12 at 20:03

Here's how I would go about it. Collect enough samples that the distribution seems stable. Then tally up the generating function for this distribution; that is, define

$$F(z) = \sum_{k=0}^\infty P(Z=k)z^k$$

Find the zeros of $$F(z)$$ by simply plotting it. Since in the limit, $$F(z)$$ should consist entirely of factors of the form $$1-p_i+p_iz$$, each of which produces a zero at $$z_i = 1-\frac{1}{p_i}$$, we can extract $$p_i = \frac{1}{1-z_i}$$.

For instance, I ran a simulation with three Bernoulli random variables. After a hundred million runs (!), I obtained the empirical distribution

$$P(Z = 0) = 0.18001324 \\ P(Z = 1) = 0.41996961 \\ P(Z = 2) = 0.32002814 \\ P(Z = 3) = 0.07998901$$

So we write

$$F(z) = 0.18001324+0.41996961z+0.32002814z^2+0.07998901z^3$$

This has three real zeros, at approximately $$-1.004, -1.435, -1.562$$, suggesting Bernoulli variables with parameters $$0.499, 0.411, 0.390$$.

In fact, the parameters were $$0.5, 0.4, 0.4$$. The twin occurrence of $$0.4$$ means that what appeared as a pair of closely spaced roots in this sampling might well appear as a close miss in another sampling. Another observation is that a hundred million samples still yielded a fair amount of error. Maybe there's another approach that's more scalable?

Here's another example with three distinct probabilities, also at a hundred million samples:

We have zeros here, empirically, at $$-0.333, -0.668, -4.000$$, implying Bernoulli variables with parameters $$0.750, 0.600, 0.200$$, and in this case, that's exactly what they were. Maybe double roots really cause a problem (unsurprisingly).

ETA: I think that if there are no doubles, the problem is computationally easier. I did a run of just a million samples, leading to a fifth-degree polynomial

$$F(z) = 0.016254+0.132128z+0.346866z^2+0.357263z^3+0.136586z^4+0.011003z^5$$

This has roots at about $$-0.2578, -0.4338, -0.9783, -1.484, -9.26$$, suggesting Bernoulli variables with parameters $$0.795, 0.697, 0.505, 0.403, 0.097$$.

As you might guess, the actual parameters were $$0.8, 0.7, 0.5, 0.4, 0.1$$, so the agreement is already pretty close after a million samples.

• Did you mean $p_i = \frac{1}{1-z_i}$? – Fred Frey Mar 12 at 21:52
• Okay, I've made a jupyter notebook that simulates this, and it seems to work okay, though it does require a lot of samples. – Fred Frey Mar 12 at 22:03
• @FredFrey: Yes, I think I must have meant that; thanks for the catch! – Brian Tung Mar 12 at 22:22
• @FredFrey: Incidentally, your error computation does not compare between ps and eps that are sorted. The second and third entries are swapped between the two. – Brian Tung Mar 12 at 22:27

For $$M=3$$, the distribution is

$$\rho_Z(0)= q_0q_1q_2,\\\rho_Z(1)= p_0q_1q_2+q_0p_1q_2+q_0q_1p_2,\\\rho_Z(2)= p_0p_1q_2+p_0q_1p_2+q_0p_1p_2,\\\rho_Z(3)= p_0p_1p_2.$$

To estimate those probabilities you can perform sufficiently many drawings to get an approximation of the $$\text{pdf}$$ of $$Z$$, then solve the non-linear system of three equations in three unknowns (the fourth equation is not independent). This generalizes to higher $$M$$.

The value of $$M$$ is simply estimated as the highest value achieved by $$Z$$, though it has a probability $$(1-p_0p_1\cdots p_M)^K$$ of not being reached after $$K$$ drawings.

• This sounds convincing to me. I will try to simulate it, and if it works I will accept the answer. – Fred Frey Mar 12 at 20:30
• I wonder if this is essentially equivalent to my approach, or if there is some deeper difference. – Brian Tung Mar 12 at 20:32