%% converted to LaTeX by EDS on 04may09

%mvp.tex:A Symbolic Computation Approach to a 
%Problem involving Multivariate Poisson distributions
%%a Plain TeX file by Eduardo Sontag Doron Zeilberger (x pages)

\documentclass[11pt]{article}

\newcommand {\EDS}[1]{\mbox{ }\medskip \\ {\bf ***** {#1} *****} \mbox{ }\medskip \\}

\usepackage{amssymb}
\usepackage{amsmath}

\topmargin -0.5in
\textheight 9.0in
\oddsidemargin 0.25in
\evensidemargin 0.25in
\textwidth 6.25in
\parskip=5pt plus 1pt minus 1pt
\parindent=0pt

\newcommand{\atopnew}[2]{\genfrac{}{}{0pt}{}{#1}{#2}}
\newcommand{\Prob}{\mbox{Pr}\,}

\title{A Symbolic Computation Approach to a Problem\\
Involving Multivariate Poisson Distributions\\
(Draft, do not distribute!)\footnote{%
Accompanied by Maple package {MVPoisson}
downloadable from 
\hfill\break
{\tt http://www.math.rutgers.edu/\~{}zeilberg/mamarim/mamarimhtml/mvp.html} }}

\author{Eduardo D.\ Sontag and Doron Zeilberger\\
Department of Mathematics\\
Rutgers University\\
Hill Center-Busch Campus\\
110 Frelinghuysen Rd., Piscataway, NJ 08854-8019, USA.\\
{\tt [sontag,zeilberg]  at math dot rutgers dot edu}}

\begin{document}

\maketitle


\section{Introduction}

Eduardo, please write it.

\EDS{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.
}




\section{The Problem}

Suppose that we have $n$ independent Poisson random variables, $X_j$ 
($j=1 \dots n$), with parameters
$\lambda_j$ respectively. In other words
\begin{equation}
\label{eq:P1}
\Prob(X_1=k_1, X_2=k_2, \dots , X_n=k_n)=
e^{-(\lambda_1 + \dots + \lambda_n)}
\frac{\lambda_1^{k_1}}{k_1!}\frac{\lambda_2^{k_2}}{k_2!}
\dots  \frac{\lambda_n^{k_n}}{k_n!} \quad .
\end{equation}

Suppose that we can't observe the $X_j$'s directly, 
but only a certain number, $m$, of linear combinations
of them:
\[
Y_i=\sum_{j=1}^{n} a_{ij} X_j \quad , \quad (i=1 , \dots , m )\quad ,
\]
where $(a_{ij})$ is a certain $m \times n$ matrix with 
non-negative coefficients.

We are interested in the following questions:

\begin{enumerate}
\item
Can one compute (fast!), for any given vector $(b_1, \dots, b_m)$ 
(possibly with large coordinates), the
probability
\[
F(b_1, \dots, b_m):= \Prob(Y_1=b_1, \dots, Y_m=b_m) \quad .
\]
\item
Can one compute (fast!), for any given vector $(b_1, \dots, b_m)$, 
(possibly with large coordinates) the
conditional expectation
\[
G_j(b_1, \dots, b_m):= 
E[X_j \,\, \bigl \vert \, \, Y_1=b_1, \dots , Y_m=b_m] \quad , \quad 
(1 \leq j \leq n) .
\]
\item
More generally, can one compute (fast!), 
the higher moments
\[
G_j^{(r)} (b_1, \dots, b_m):= 
E[X_j^r \,\, \bigl \vert \,\, Y_1=b_1, \dots , Y_m=b_m]
\quad , \quad (r \geq 2) \quad ,
\]
that would immediately allow us to compute the moments about the mean.
Can we compute (fast!) mixed moments, in particular the covariances?
\end{enumerate}

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

\EDS{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?}

\section{The generating function}

Fix a matrix $A=(a_{ij})$ ($1 \leq i \leq m$,  $1 \leq j \leq n$),
once and for all. 
Let
\[
F_0(b_1, \dots, b_m)=
%e^{-(\lambda_1 + \dots + \lambda_n)}
\sum_{
\atopnew{k_1, \dots, k_n \geq 0}%
{a_{11} k_1 + \dots + a_{1n} k_n=b_1 , \dots , 
a_{m1} k_1 + \dots + a_{mn} k_n=b_m}
}
\frac{\lambda_1^{k_1}}{k_1!}\frac{\lambda_2^{k_2}}{k_2!} \dots  
\frac{\lambda_n^{k_n}}{k_n!}
\]
(value is zero if the sum is empty), so that
\[
F(b_1, \dots, b_m)=
e^{-(\lambda_1 + \dots + \lambda_n)}F_0(b_1, \dots, b_m) \,.
\]
Let 
%$f(z_1, \dots, z_m)$ 
$f_0$
be the (multivariable) generating function of 
$F_0$,
%$F(b_1, \dots, b_m)$,
in other words
\[
f_0(z_1, \dots, z_m)=\sum_{b_1 \geq 0, \dots, b_m \geq 0} 
F_0(b_1, \dots, b_m) z_1^{b_1} \cdots z_m^{b_m} \quad.
\]

%%%THIS PART IS NOT NEEDED. We should phrase everying in terms of F_0 and f_0
%Thus,
%\[
%f(z_1, \dots, z_m)=\sum_{b_1 \geq 0, \dots, b_m \geq 0} 
%F(b_1, \dots, b_m) z_1^{b_1} \cdots z_m^{b_m}
%\]
%is the generating function of $F$.
%%%End of Deleted part

Our quantity of interest, $F_0(b_1, \dots, b_m)$,
is the coefficient of $z_1^{b_1} \cdots z_m^{b_m}$ in the multivariable
Taylor expansion about the origin of $f_0(z_1, \dots, z_m)$.

We have:
\[
%e^{\lambda_1 + \dots + \lambda_n}
f_0(z_1, \dots, z_m)
%=\sum_{b_1 \geq 0, \dots, b_m \geq 0} F_0(b_1, \dots, b_m) z_1^{b_1} \cdots z_m^{b_m} =
=\sum_{b_1 \geq 0, \dots, b_m \geq 0} 
\left (
\sum_{
\atopnew{k_1, \dots, k_n \geq 0}%
{a_{11} k_1 + \dots + a_{1n} k_n=b_1 , \dots , 
a_{m1} k_1 + \dots + a_{mn} k_n=b_m}
}
\frac{\lambda_1^{k_1}}{k_1!}\frac{\lambda_2^{k_2}}{k_2!} \dots  
\frac{\lambda_n^{k_n}}{k_n!}
\right )
z_1^{b_1} \cdots z_m^{b_m} \quad .
\]
By changing the order of summation, this equals
\newcommand{\beqn}{\begin{eqnarray*}}
\newcommand{\eeqn}{\end{eqnarray*}}
\beqn
&& \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!
\sum_{k_1 \geq 0, \dots, k_n \geq 0}
\frac{\lambda_1^{k_1}}{k_1!}\frac{\lambda_2^{k_2}}{k_2!} \dots  \frac{\lambda_n^{k_n}}{k_n!}
z_1^{a_{11}k_1+ \dots + a_{1n} k_n} \cdots z_m^{a_{m1} k_1 + \dots + a_{mn} k_n}
\\
&=&\sum_{k_1 \geq 0, \dots, k_n \geq 0} 
\frac{(\lambda_1 z_1^{a_{11}} z_2^{a_{21}} \cdots z_m^{a_{m1}})^{k_1}}{k_1!}
\cdots
\frac{(\lambda_n z_1^{a_{1n}} z_2^{a_{2n}} \cdots z_m^{a_{mn}})^{k_n}}{k_n!}
\\
&=&
\left ( \sum_{k_1 \geq 0}  \frac{(\lambda_1 z_1^{a_{11}} z_2^{a_{21}} \cdots z_m^{a_{m1}})^{k_1}}{k_1!} \right )
\cdots
\left ( \sum_{k_n \geq 0}  \frac{(\lambda_n z_1^{a_{1n}} z_2^{a_{2n}} \cdots z_m^{a_{mn}})^{k_n}}{k_n!} \right )
\\
&=&\exp ( \lambda_1 z_1^{a_{11}} z_2^{a_{21}} \cdots z_m^{a_{m1}} )
\cdots
\exp ( \lambda_n z_1^{a_{1n}} z_2^{a_{2n}} \cdots z_m^{a_{mn}} )
\\
&=&\exp \left ( \lambda_1 z_1^{a_{11}} z_2^{a_{21}} \cdots z_m^{a_{m1}} + \cdots +
\lambda_n z_1^{a_{1n}} z_2^{a_{2n}} \cdots z_m^{a_{mn}} \right ) \quad .
\eeqn

We have just derived

{\bf Theorem 1}: 
\[
f_0(z_1, \dots, z_m)= 
%%%Delete this line: e^{-(\lambda_1 + \dots + \lambda_n)}
\exp \left ( 
\sum_{j=1}^{n} \lambda_j \prod_{i=1}^{m} z_i^{a_{ij}} 
\right )
\]

The conditional probability
\[
\Pr\left(X_1=k_1, X_2=k_2, \dots , X_n=k_n
        \,\, \bigl \vert \, \, Y_1=b_1, \dots , Y_m=b_m\right)
\]
is the same as the expression in (\ref{eq:P1})
%% label the first display in the current writeup as {eq:P1}
divided by $F(b)$, provided that $\sum_{j=1}^{n} a_{ij} k_j=b_i$ for all $i$,
and is zero otherwise.
Recall that the $r$th 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)}\, \bigl \vert \, Y]$.
By definition,
$E[X_j^{(r)}\, \bigl \vert \, Y]$
%%%Delete this line equals $e^{-(\lambda_1 + \dots + \lambda_n)}$ multiplied by
the following expression and divided by $F_0(b)$: %%%Changed F(b) to F_0(b) becuause we chucked the e^{-( ....)
\begin{equation}
\label{eq:moment1}
\sum_{
\atopnew{k_1, \dots, k_n \geq 0}
{a_{11} k_1 + \dots + a_{1n} k_n=b_1 , \dots ,
a_{m1} k_1 + \dots + a_{mn} k_n=b_m}}
k_j(k_j-1)\dots (k_j-r+1)
\frac{\lambda_1^{k_1}}{k_1!}\frac{\lambda_2^{k_2}}{k_2!} \dots
\frac{\lambda_n^{k_n}}{k_n!} \quad .
% will label this as eq:moment1
\end{equation}
Now, expression (\ref{eq:moment1}) is the same as
the result of applying the operator
$\lambda_j^r (\frac{\partial}{\partial \lambda_j})^r$ to
$F_0(b_1, \dots, b_m)$
\emph{when viewing the $\lambda$'s as variables and not as constants}.
On the other hand,
\[
\lambda_j^r (\frac{\partial}{\partial \lambda_j})^r f_0(z_1,\dots,z_m)
=
\sum_{b_1 \geq 0, \dots, b_m \geq 0}
\lambda_j^r (\frac{\partial}{\partial \lambda_j})^r
F_0(b_1, \dots, b_m) z_1^{b_1} \cdots z_m^{b_m}
\]
and therefore expression (\ref{eq:moment1}) is the same as the
coefficient of $z_1^{b_1} \cdots z_m^{b_m}$ in
$\lambda_j^r (\frac{\partial}{\partial \lambda_j})^r f_0(z_1,\dots,z_m)$
Since, as formal power series, we have the representation in Theorem 1.
%% to replace theorem "1" by soft label
we conclude that expression (\ref{eq:moment1}) is the same as
the coefficient of
$z_1^{b_1} \cdots z_m^{b_m}$ in
$(\prod_{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}, \dots, b_m-ra_{mj})$
when all $b_i-ra_{ij}\geq 0$ and zero otherwise.
In conclusion, $E[X_j^{(r)}\, \bigl \vert \, Y]$
equals 
$F_0(b_1-ra_{1j},b_2- ra_{2j}, \dots, b_m-ra_{mj})$
divided by $F_0(b)$.
We have proved:

{\bf Theorem 2}: The conditional factorial moments
$E[X_j^{(r)}\, \bigl \vert \, Y]$ are given in terms of the
$F_0(b_1, \dots , b_m)$ by
\[
\lambda_j^r \, \cdot \,
\frac{F_0(b_1-ra_{1j},b_2- ra_{2j}, \dots, b_m-ra_{mj}) }%
{F_0(b_1, \dots, b_m)}
\]
when all $b_i-ra_{ij}\geq 0$ and zero otherwise.

So everything depends on a fast computation of the 
coefficients $F_0(b_1, \dots , b_m)$,  of $f_0(z_1, \dots , z_m)$.


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

\section{Recurrences}

From now on, let's assume that the entries of $A$, $(a_{ij})$,
are non-negative {\bf integers}. In that case, we can write
\EDS{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, \dots, z_m) ) \quad,
\]
where $Q(z_1, \dots , z_m)$ is the {\bf polynomial}
\[
Q(z_1, \dots, z_m) :=
\sum_{j=1}^{n} \lambda_j \prod_{i=1}^{m} z_i^{a_{ij}} 
\quad .
\]
By Cauchy's theorem, we can express $F(b_1, \dots, b_m)$ as
a {\bf multi-contour integral}:
\[
F(b_1, \dots, b_m) =
\left ( \frac{1}{2 \pi i}  \right )^m
\int_{|z_1|=c} \dots \int_{|z_m|=c}
\frac {\exp( Q(z_1, \dots, z_m) )} {z_1^{b_1+1} \cdots z_m^{b_m+1} }
\, dz_1 \, \dots \, dz_m
\quad .
\]
By the celebrated {\bf Wilf-Zeilberger} theory (\cite{[WZ]}),
$F(b_1, \dots, b_m)$ satisfies pure {\bf 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^{(i)}_r(b_1, \dots, b_m)$ ($0 \leq r \leq R_i$)
such that the following
holds, for \emph{all} $(b_1, \dots , b_m)$:
\[
\sum_{r=0}^{R_i} P^{(i)}_r(b_1, \dots, b_m)
F(b_1, \dots, b_{i-1}, b_i+r, b_{i+1} , \dots , b_m) = 0 \quad .
\]
Once these recurrences are known, one can compute
$F(b_1, \dots, b_m)$ in 
time linear in $b_1 + \dots + b_m$ and 
with constant memory allocation (one only needs to
remember, at each stage, a constant number of values).
\EDS{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 {\bf Apagodu-Zeilberger}\cite{[ApZ]}
multi-variable extension of the Almkvist-Zeilberger\cite{[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.  

%{\vfill \eject}

\section{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 \emph{single sum}. Indeed, let

%$c_{01}$ be the some of $g_j$ for which $a_{1,j}=0$, $a_{2,j}=1$,
%
%$c_{10}$ be the some of $g_j$ for which $a_{1,j}=1$, $a_{2,j}=0$,
%
%$c_{11}$ be the some of $g_j$ for which $a_{1,j}=1$, $a_{2,j}=1$, then we have

\begin{itemize}
\item
$c_{01}$ be the sum of the $\lambda_j$'s for which $a_{1,j}=0$, $a_{2,j}=1$,
\item
$c_{01}$ be the sum of the $\lambda_j$'s  for which $a_{1,j}=1$, $a_{2,j}=0$,
\item
$c_{01}$ be the sum of the $\lambda_j$'s for which $a_{1,j}=1$, $a_{2,j}=1$.
\end{itemize}
Then, we have
\[
Q(z)=c_{01} z_1+c_{10}z_2+c_{11}z_1z_2 ,
\]
and so
%\EDS{must watch for missing exponential, as earlier..., DZ: yes, f should be replaced by f_0}
\beqn
f_0(z_1,z_2) &=& e^{Q(z)} \;=\; \sum_{k=0}^{\infty} \frac{Q(z)^k}{k!}
\\
&=&
\sum_{\alpha \geq 0, \beta \geq 0, \gamma \geq 0} 
 \frac{ (c_{01} z_1)^{\alpha} (c_{10} z_2)^{\beta} (c_{11} z_1 z_2)^{\gamma}}%
      {\alpha! \beta! \gamma!}
\\
&=&\sum_{\alpha \geq 0, \beta \geq 0, \gamma \geq 0} 
 \frac{ c_{01}^{\alpha} c_{10}^{\beta} c_{11}^{\gamma} z_1^{\alpha+\gamma} z_2^{\beta+\gamma} }%
  {\alpha! \beta! \gamma!} \quad .
\eeqn
To get $F_0(b_1,b_2)$, we must extract the coefficient of $z_1^{b_1} z_2^{b_2}$
 which entails
$\alpha=b_1-\gamma,\beta=b_2-\gamma$, and we have the single-sum binomial coefficient (hypergeometric) sum
(replacing $\gamma$ by $k$)
\[
F_0(b_1,b_2)=\sum_{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)!} \quad .
\]
Using the {\bf Zeilberger Algorithm} (\cite{[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_1b_2)F_0(b_1+2,b_2)
+
\]
\[
(-c_{11}b_1-c_{11}+b_2c_{11}+b_2c_{10}c_{01}-2b_1c_{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 \quad .
\]


\section{The Maple package MVPoisson}

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

{\tt http://www.math.rutgers.edu/\~{}zeilberg/mamarim/mamarimhtml/mvp.html} \quad ,

where one can also find sample input and output.

\EDS{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
$\Omega(n) = c n^p$?}

\EDS{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?}

\begin{thebibliography}{10}

\bibitem[AlZ]{[AlZ]}
G. Almkvist and D. Zeilberger,
\emph{The method of differentiating under the
integral sign}, J. Symbolic Computation {\bf 10}(1990), 571-591.
 
\bibitem[ApZ]{[ApZ]}
M. Apagodu  and D. Zeilberger,
\emph{Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms and the
Sharpening of Wilf-Zeilberger Theory },
Adv. Appl. Math. {\bf 37} (2006)(Special Regev issue), 139-152

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

\bibitem[WZ]{[WZ]}
H.S. Wilf and D. Zeilberger,
\emph{An algorithmic proof theory for hypergeometric
(ordinary and ``q") multisum/integral identities}, Invent. Math. 
{\bf 108} (1992), 575-633.

\bibitem[Z]{[Z]}
D. Zeilberger,
\emph{The method of creative telescoping}, J. Symbolic Computat.
{\bf 11}, 195-204 (1991).

\end{thebibliography}

\end{document}
