A Symbolic Computation Approach to a Problem Involving Multivariate Poisson Distributions (Draft, do not distribute!)^{∗}

Eduardo D. Sontag and Doron Zeilberger Department of Mathematics Rutgers University Hill Center-Busch Campus 110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.

              [sontag,zeilberg] at math dot rutgers dot edu

May 5, 2009

1 Introduction

Eduardo, please write it.

***** important to decide: do we really start by emphasizing chemistry, or do we first state the purely mathematical problem (including Prussian officers and email messages. well, no Prussian officers, but we can make a good story with text messages, emails, and so on)? Then the chemistry would be a self-contained second part. DZ’s reply: Yes, good idea! Keep the Chemistry to the second part, but mention it in the still-to-be-written abstart. *****

2 The Problem

Suppose that we have n independent Poisson random variables, X_{j} (j = 1 … n), with parameters λ_{j} respectively. In other words

Pr (X_{1} = k_{1}, X_{2} = k_{2}, …, X_{n} = k_{n}) = e^{−(λ_{1}+···+λ_{n})} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} .  (1)

∗Accompanied by Maple package MVPoisson downloadable from http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/mvp.html

Suppose that we can’t observe the X_{j}’s directly, but only a certain number, m, of linear combinations of them:

Y_{i} = ∑_{j=1}^{n} a_{ij}X_{j} , (i = 1, …, m) ,

where (a_{ij}) is a certain m × n matrix with non-negative coefficients.

We are interested in the following questions:

1. Can one compute (fast!), for any given vector (b_{1}, …, b_{m}) (possibly with large coordinates), the probability

F (b_{1}, …, b_{m}) := Pr (Y_{1} = b_{1}, …, Y_{m} = b_{m}) .

2. Can one compute (fast!), for any given vector (b_{1}, …, b_{m}), (possibly with large coordinates) the conditional expectation

G_{j}(b_{1}, …, b_{m}) := E[X_{j} | Y_{1} = b_{1}, …, Y_{m} = b_{m}] , (1 ≤ j ≤ n).

3. More generally, can one compute (fast!), the higher moments

G_{j}^{(r)}(b_{1}, …, b_{m}) := E[X_{j}^{r} | Y_{1} = b_{1}, …, Y_{m} = b_{m}] , (r ≥ 2) ,

that would immediately allow us to compute the moments about the mean. Can we compute (fast!) mixed moments, in particular the covariances?

***** EDS will fill-in: As trivial example, consider the case in which there is one column... multinomial *****

***** actually, we are also VERY interested in covariances, so we need JOINT moments too. Doron: could you write that up easily, or should we say something like “one can also easily find them” for now, if we want to be done soon? *****

3 The generating function

Fix a matrix A = (a_{ij}) (1 ≤ i ≤ m, 1 ≤ j ≤ n), once and for all. Let

F_{0}(b_{1}, …, b_{m}) = ∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+···+a_{1n}kn=b_{1},...,a_{m1}k_{1}+···+amnkn=bm}} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!}

(value is zero if the sum is empty), so that

F (b_{1}, …, b_{m}) = e^{−(λ_{1}+···+λ_{n})}F_{0}(b_{1}, …, b_{m}) .

Let f_{0} be the (multivariable) generating function of F_{0}, in other words

f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} F_{0}(b_{1}, …, b_{m})z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}} .

Our quantity of interest, F_{0}(b_{1}, …, b_{m}), is the coefficient of z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}} in the multivariable Taylor expansion about the origin of f_{0}(z_{1}, …, z_{m}).

We have:

f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} ( ∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+···+a_{1n}kn=b_{1},...,a_{m1}k_{1}+···+amnkn=bm}} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} ) z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}} .

By changing the order of summation, this equals

∑_{k_{1}≥0,...,k_{n}≥0} \frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} z_{1}^{a_{11}k_{1}+···+a_{1n}k_{n}} ⋯ z_{m}^{a_{m1}k_{1}+···+a_{mn}k_{n}}

= ∑_{k_{1}≥0,...,k_{n}≥0} \frac{(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} ⋯ z_{m}^{a_{m1}})^{k_{1}}}{k_{1}!} ⋯ \frac{(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} ⋯ z_{m}^{a_{mn}})^{k_{n}}}{k_{n}!}

= ( ∑_{k_{1}≥0} \frac{(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} ⋯ z_{m}^{a_{m1}})^{k_{1}}}{k_{1}!} ) ⋯ ( ∑_{k_{n}≥0} \frac{(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} ⋯ z_{m}^{a_{mn}})^{k_{n}}}{k_{n}!} )
= exp(λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} ⋯ z_{m}^{a_{m1}}) ⋯ exp(λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} ⋯ z_{m}^{a_{mn}})
= exp (λ_{1}z_{1}^{a_{11}}z_{2}^{a_{21}} ⋯ z_{m}^{a_{m1}} + ⋯ + λ_{n}z_{1}^{a_{1n}}z_{2}^{a_{2n}} ⋯ z_{m}^{a_{mn}}) .

We have just derived

Theorem 1:

f_{0}(z_{1}, …, z_{m}) = exp (∑_{j=1}^{n} λ_{j} ∏_{i=1}^{m} z_{i}^{a_{ij}} )

The conditional probability

Pr (X_{1} = k_{1}, X_{2} = k_{2}, …, X_{n} = k_{n} | Y_{1} = b_{1}, …, Y_{m} = b_{m})

is the same as the expression in (1) divided by F (b), provided that ∑_{j=1}^{n} a_{ij}k_{j} = b_{i} for all i, and is zero otherwise. Recall that the rth factorial moment of a random variable W, E[W^{(r)}], is, by definition, the expectation of W!/(W − r)!. We are interested in the conditional factorial moments of X_{j} given Y = b, which we will denote as E[X_{j}^{(r)} | Y]. By definition, E[X_{j}^{(r)} | Y] the following expression and divided by F_{0}(b):

∑_{{k_{1},...,kn≥0 \atop a_{11}k_{1}+···+a_{1n}kn=b_{1},...,a_{m1}k_{1}+···+amnkn=bm}} k_{j}(k_{j} − 1) … (k_{j} − r + 1)\frac{λ_{1}^{k_{1}}}{k_{1}!} \frac{λ_{2}^{k_{2}}}{k_{2}!} … \frac{λ_{n}^{k_{n}}}{k_{n}!} .  (2)

Now, expression (2) is the same as the result of applying the operator λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r} to F_{0}(b_{1}, …, b_{m}) when viewing the λ’s as variables and not as constants. On the other hand,

λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}f_{0}(z_{1}, …, z_{m}) = ∑_{b_{1}≥0,...,b_{m}≥0} λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}F_{0}(b_{1}, …, b_{m})z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}}

and therefore expression (2) is the same as the coefficient of z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}} in λ_{j}^{r}(\frac{∂}{∂λ_{j}})^{r}f_{0}(z_{1}, …, z_{m}) Since, as formal power series, we have the representation in Theorem 1. we conclude that expression (2) is the same as the coefficient of z_{1}^{b_{1}} ⋯ z_{m}^{b_{m}} in (∏_{i=1}^{m} z_{i}^{a_{ij}})^{r}f_{0}(z), which is the same as F (b_{1} −ra_{1j}, b_{2} −ra_{2j}, …, b_{m} −ra_{mj}) when all b_{i} −ra_{ij} ≥ 0 and zero otherwise. In conclusion, E[X_{j}^{(r)} | Y] equals F_{0}(b_{1} − ra_{1j}, b_{2} − ra_{2j}, …, b_{m} − ra_{mj}) divided by F_{0}(b). We have proved:

Theorem 2: The conditional factorial moments E[X_{j}^{(r)} | Y] are given in terms of the F_{0}(b_{1}, …, b_{m}) by

λ_{j}^{r} · \frac{F_{0}(b_{1} − ra_{1j}, b_{2} − ra_{2j}, …, b_{m} − ra_{mj})}{F_{0}(b_{1}, …, b_{m})}

when all b_{i} − ra_{ij} ≥ 0 and zero otherwise.

So everything depends on a fast computation of the coefficients F_{0}(b_{1}, …, b_{m}), of f_{0}(z_{1}, …, z_{m}).

By taking mixed partial derivatives, we can easily derive analogous expressions for mixed moments, in particular, the covariances.

4 Recurrences

From now on, let’s assume that the entries of A, (a_{ij}), are non-negative integers. In that case, we can write

***** next line is not quite right unless the exponentials are taken into account in the definition of Q. Or, equivalently, we can compute f_{0}. DZ: Indeed we should stick to f_{0}. *****

f_{0}(z) = exp(Q(z_{1}, …, z_{m})) ,

where Q(z_{1}, …, z_{m}) is the polynomial

Q(z_{1}, …, z_{m}) := ∑_{j=1}^{n} λ_{j} ∏_{i=1}^{m} z_{i}^{a_{ij}} .

By Cauchy’s theorem, we can express F (b_{1}, …, b_{m}) as a multi-contour integral:

F (b_{1}, …, b_{m}) = (\frac{1}{2πi})^{m} ∫_{|z_{1}|=c} … ∫_{|z_{m}|=c} \frac{exp(Q(z_{1}, …, z_{m}))}{z_{1}^{b_{1}+1} ⋯ z_{m}^{b_{m}+1}} dz_{1} … dz_{m} .

By the celebrated Wilf-Zeilberger theory ([WZ]), F (b_{1}, …, b_{m}) satisfies pure linear recurrences with polynomial coefficients in each of its arguments. This means that for each i between 1 and m, there exists a positive integer R_{i} (the order) and polynomials P_{r}^{(i)}(b_{1}, …, b_{m}) (0 ≤ r ≤ R_{i}) such that the following holds, for all (b_{1}, …, b_{m}):

∑_{r=0}^{R_{i}} P_{r}^{(i)}(b_{1}, …, b_{m})F (b_{1}, …, b_{i−1}, b_{i} + r, b_{i+1}, …, b_{m}) = 0 .

Once these recurrences are known, one can compute F (b_{1}, …, b_{m}) in time linear in b_{1} +⋯+b_{m} and with constant memory allocation (one only needs to remember, at each stage, a constant number of values).

***** add a warning about zeros of coefs of recursion? DZ: OK, here goes *****

In rare cases, the leading term of the recurrence would vanish, in which case, we would encounter a (discrete) “singularity”, and would not be able to go on, since we would have to divide by 0, but in that case one can show that there is an alternative route, using another order of applying the recurrences.

The Apagodu-Zeilberger[ApZ] multi-variable extension of the Almkvist-Zeilberger[AlZ] algorithm can find such recurrences explicitly. Unfortunately, for matrices A with more than three rows, the time taken to find such recurrences is prohibitive, but many matrices of interest have two or three rows.

5 Two-Rowed matrices

If the matrix A only has two rows, and the entries are only {0, 1}, then one can express F_{(}b_{1}, b_{2}) as a single sum. Indeed, let

• c_{01} be the sum of the λ_{j}’s for which a_{1,j} = 0, a_{2,j} = 1,

• c_{01} be the sum of the λ_{j}’s for which a_{1,j} = 1, a_{2,j} = 0,

• c_{01} be the sum of the λ_{j}’s for which a_{1,j} = 1, a_{2,j} = 1.

Then, we have

Q(z) = c_{01}z_{1} + c_{10}z_{2} + c_{11}z_{1}z_{2},

and so

f_{0}(z_{1}, z_{2}) = e^{Q(z)} = ∑_{k=0}^{∞} \frac{Q(z)^{k}}{k!}

= ∑_{α≥0,β≥0,γ≥0} \frac{(c_{01}z_{1})^{α}(c_{10}z_{2})^{β}(c_{11}z_{1}z_{2})^{γ}}{α!β!γ!}

= ∑_{α≥0,β≥0,γ≥0} \frac{c_{01}^{α}c_{10}^{β}c_{11}^{γ}z_{1}^{α+γ}z_{2}^{β+γ}}{α!β!γ!} .

To get F_{0}(b_{1}, b_{2}), we must extract the coefficient of z_{1}^{b_{1}}z_{2}^{b_{2}} which entails α = b_{1} − γ, β = b_{2} − γ, and we have the single-sum binomial coefficient (hypergeometric) sum (replacing γ by k)

F_{0}(b_{1}, b_{2}) = ∑_{k=0}^{min(b_{1},b_{2})} \frac{c_{11}^{k}c_{01}^{b_{1}−k}c_{10}^{b_{2}−k}}{k!(b_{1} − k)!(b_{2} − k)!} .

Using the Zeilberger Algorithm ([Z, PWZ]) , we get the following linear recurrence:

(c_{10}b_{1}^{2} + 4c_{10} − 2c_{10}b_{2} + 4c_{10}b_{1} − c_{10}b_{1}b_{2})F_{0}(b_{1} + 2, b_{2})+

(−c_{11}b_{1}−c_{11}+b_{2}c_{11}+b_{2}c_{10}c_{01}−2b_{1}c_{10}c_{01}−3c_{01}c_{10})F_{0}(b_{1}+1, b_{2})+(c_{11}c_{10}+c_{01}c_{10}^{2})F_{0}(b_{1}, b_{2}) = 0 .

6 The Maple package MVPoisson

All this is implemented in the Maple package MVPoisson accompanying this article. It is available from the webpage of this article

http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/mvp.html ,

where one can also find sample input and output.

***** paper: I looked at the sample outputs on the web, and saw statements such as:

“This seems to be asymptotically etc etc...”

Can you please include a sentence explaining how this fit is done by the software? Or did you do this by hand? Is there any reason to expect Ω(n) = cn^{p}? *****

***** I’d like to download all at once all the files linked from

http://www.math.rutgers.edu/˜ zeilberg/mamarim/mamarimhtml/mvp.html, but for some strange reason my “wget” is failing to work recursively (you must have a funny directory structure). Is it possible to get a tar file? *****

References

[AlZ] G. Almkvist and D. Zeilberger, The method of differentiating under the integral sign, J. Symbolic Computation 10(1990), 571-591.

[ApZ] M. Apagodu and D. Zeilberger, Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the Sharpening of Wilf-Zeilberger Theory , Adv. Appl. Math. 37 (2006)(Special Regev issue), 139-152

[PWZ] M. Petkovsek, H.S. Wilf and D. Zeilberger, A=B, AK Peters, Wellesley, (1996). [available on-line from the authors’ websites.]

[WZ] H.S. Wilf and D. Zeilberger, An algorithmic proof theory for hypergeometric (ordinary and “q”) multisum/integral identities, Invent. Math. 108 (1992), 575-633.

[Z] D. Zeilberger, The method of creative telescoping, J. Symbolic Computat. 11, 195-204 (1991).