Pushing for the Extreme: Estimation of Poisson Distribution from Low Count Unreplicated Data—How Close Can We Get?

Abstract

Studies of learning algorithms typically concentrate on situations where potentially ever growing training sample is available. Yet, there can be situations (e.g., detection of differentially expressed genes on unreplicated data or estimation of time delay in non-stationary gravitationally lensed photon streams) where only extremely small samples can be used in order to perform an inference. On unreplicated data, the inference has to be performed on the smallest sample possible—sample of size 1. We study whether anything useful can be learnt in such extreme situations by concentrating on a Bayesian approach that can account for possible prior information on expected counts. We perform a detailed information theoretic study of such Bayesian estimation and quantify the effect of Bayesian averaging on its first two moments. Finally, to analyze potential benefits of the Bayesian approach, we also consider Maximum Likelihood (ML) estimation as a baseline approach. We show both theoretically and empirically that the Bayesian model averaging can be potentially beneficial.

1. Introduction

Studies in (computational) learning theory mostly tend to concentrate on situations where potentially ever increasing number of training examples is available. While such results can lead to deep insights into the workings of learning algorithms, e.g., linking together characteristics of the data generating distributions, learning machines and sample sizes, there can be situations where, by very nature of the problem, only extremely small samples are available. In such situations it is of utmost importance to theoretically analyze exactly what and under what circumstances can be learnt. One example of such a scenario in count data is detection of differentially expressed genes, where even subtle changes in gene expression levels can be indicators of biologically crucial processes [1]. When replicas are costly to obtain one can attempt to use the limited data at one’s disposal to make the relevant inferences, as for example in the Audic and Claverie approach [2,3,4,5,6]. Another situation where available count data can be extremely sparse is estimation of time delay in non-stationary gravitationally lensed photon streams. When the scale of variability of the source is of order, say, of more than tens of days and observation gaps are not too long, one can resolve the time delay between lensed images of the same source by working directly with daily measurements of fluxes in the radio, optical or X-ray range [7,8,9,10]. However, when the variability scale is of the order of hours, one must turn to photon streams in the lensed images. One possibility of time delay detection in such cases is through comparing counts in relatively short and time-shifted moving time windows in the lensed photon streams.

In this paper we theoretically study what happens in the extreme situation of unreplicated data when the inference has to be performed on the smallest sample possible—sample of size 1. We consider a model-based Bayesian approach that averages over possible Poisson models with weighting determined by the posterior over the models, given the single observation. In fact, such a Bayesian approach has been considered in the bioinformatics literature under the assumption of flat improper prior over the Poisson rate parameter [2,3,4,5,6]. One can, of course, be excused for being highly sceptical about the relevance of such inferences, yet the methodology has apparently been used in a number of successful studies. In an attempt to build theoretical foundations behind such inference schemes, we proved a rather surprising result [11]: The expected Kullback–Leibler divergence from the true unknown Poisson distribution to its model learnt from a single realization never exceeds 1/2 bit.

Even though the field of bioinformatics is moving fast and better procedures for detection of differentially expressed genes have been introduced (e.g., not relying on the Poisson assumption, specifically taking into account potential dependencies among the genes, etc.), the primary focus of this study is different. Irrespective of the application domain, we theoretically investigate how reliably can a model for count data be build from a single count observation, under the assumption of a Poisson source. There are two issues that need careful consideration:

• Equal a-priori weighting (flat prior) over possible (unknown) Poisson sources is unrealistic. Typical values of observed counts are usually bounded by the nature of the problem (e.g., gene magnification setting used in the experiments or time window on the photon streams). One may have a good initial (a-priori) guess as to what ranges of typical observed counts might be reasonably expected. In particular, we are interested in the low count regimes. In such cases, it is desirable to incorporate such prior knowledge into the inference mechanism. In this study, we do this in the Bayesian framework through prior distribution over the expected counts.
• To understand potential benefits of the proposed learning/inference method (in our case Bayesian approach), it is important to compare it with a simple straightforward baseline (here maximum likelihood estimation). We contrast the expected Kullback–Leibler divergences from the true unknown Poisson distribution to its Bayesian and maximum likelihood estimates, inferred from a single realization.

The paper has the following organization. In Section 2 we introduce the maximum likelihood and Bayesian (with flat prior over mean rates) approaches to inferring predictive distribution over counts based on a single count observation. We also briefly review past work on information theoretic properties of the two approaches. Section 3 contains derivation of a more general Bayesian approach with gamma prior on the mean count parameter. In Section 4 we calculate the first two central moments of our generalized model. This enables us to better understand the influence of the prior on the inferred model and highlight the differences with the previous approach using the flat (improper) prior. In Section 5 we perform an information theoretic study of learning capabilities of the generalized model. Empirical investigations are presented in Section 6 and the main findings are discussed and summarized in Section 7.

2. Single Count Data—Bayesian and Maximum Likelihood Approaches

In this section we will briefly review the original Audic–Claverie [2] and maximum likelihood approaches outside the bioinformatics context.

2.1. Bayesian Averaging in the Audic–Claverie Approach

Let x be an observed count in an experiment. When repeating the experiment, possibly under different conditions, we observe a (possibly different) count y. The quantity of interest is the probability of observing y given that we already observed x, not knowing the identity of the generating Poisson source

$$P(X=x|\lambda )=e^{-\lambda }\frac{{\lambda }^x}{x!}$$

(1)

where $\lambda \ge 0$ is the (unknown) parameter representing the mean count value.

Under the null hypothesis (not differentially expressed genes), both counts x and y come from the same underlying Poisson distribution $P(·|\lambda )$. The key instrument in the Audic–Claverie approach is a distribution $P_{AC}\left(y\right|x)$ over counts y informed by the observed count x, under the null hypothesis. $P_{AC}\left(y\right|x)$ is obtained by Bayesian averaging (infinite mixture) of all possible Poisson distributions $P\left(y\right|{\lambda }^{\prime })$ with mixing proportions equal to the posteriors $p\left({\lambda }^{\prime }\right|x)$ under the flat prior over λ. Formally, the probability of count y, given the observed count x from the same (unknown) Poisson distribution is:

$$\begin{array}{ccc}\hfill P\left(y\right|x)& =& {\int }_0^{\infty }p(y,\lambda |x)\text{d}\lambda \hfill \\ & =& {\int }_0^{\infty }P\left(y\right|\lambda ,x)p\left(\lambda \right|x)\text{d}\lambda \hfill \\ & =& {\int }_0^{\infty }P\left(y\right|\lambda )\frac{P\left(x\right|\lambda \left)p\right(\lambda )}{{\int }_0^{\infty }P\left(x\right|{\lambda }^{\prime })p\left({\lambda }^{\prime }\right)\text{d}{\lambda }^{\prime }}\text{d}\lambda \hfill \end{array}$$

(2)

Imposing the flat (improper) prior $p\left(\lambda \right)$ over the Poisson parameter λ results in

$$P_{AC}\left(y\right|x)=\frac{1}{y!}\frac{{\int }_0^{\infty }{e}^{-2\lambda }{\lambda }^{x+y}\text{d}\lambda }{{\int }_0^{\infty }{e}^{-\lambda }{\lambda }^x\text{d}\lambda }$$

Since Gamma distribution parameterized by $a,b>0$ takes the form

$$Gamma\left(\lambda \right|a,b)=\frac{1}{\Gamma \left(a\right)}{b}^a{\lambda }^{a-1}{e}^{-b\lambda }$$

where $\Gamma \left(a\right)={\int }_0^{\infty }{u}^{a-1}{e}^{-u}\text{d}u$ is the Gamma function, we have

$$P_{AC}\left(y\right|x)=\frac{1}{y!2^{x+y+1}}\frac{\Gamma (x+y+1)}{\Gamma (x+1)}$$

(3)

which, since x and y are integers (i.e., $\Gamma \left(x\right)=(x-1)!$), can be rewritten as

$$\begin{array}{ccc}\hfill P_{AC}\left(y\right|x)& =& \frac{1}{2^{x+y+1}}\frac{(x+y)!}{x!y!}\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill \phantom{P_{AC}\left(y\right|x)}& =& \frac{1}{2^{x+y+1}}\left(\genfrac{}{}{0pt}{}{x+y}{x}\right)\hfill \end{array}$$

(5)

$P_{AC}(·|x)$ can then be used, e.g., for principled inferences, construction of confidence intervals or statistical testing.

2.2. Information Theory of $P_{AC}\left(y\right|x)$

Consider a “true” underlying Poisson distribution $P\left(y\right|\lambda )$ (1) over possible counts $y\ge 0$. We first use $P(·|\lambda )$ to generate a count x and then employ $P_{AC}\left(y\right|x)$ (5) as a model distribution over y, given the already observed count x. We ask: If we repeated the process above, how different, in terms of Kullback–Leibler divergence, are on average the two distributions over y? One would naturally hope that $P_{AC}\left(y\right|x)$ is sufficiently representative of the true unknown distribution $P\left(y\right|\lambda )$.

In [11] we proved that, given an underlying Poisson distribution $P\left(x\right|\lambda )$, if we repeatedly generated a “representative” count x from $P\left(x\right|\lambda )$, the average divergence $\mathcal{E}\left(\lambda \right)$ of $P_{AC}\left(y\right|x)$ from the truth $P\left(y\right|\lambda )$ would never exceed 1/2 bit.

Theorem 1 [11] Consider an underlying Poisson distribution $P(·|\lambda )$ parameterized by some $\lambda >0$. Then

$$\mathcal{E}\left(\lambda \right)=E_{P\left(x\right|\lambda )}\left[D_{KL}[P\left(y\right|\lambda )\parallel P_{AC}\left(y\right|x)]\right]=\frac{1}{2}log2+O\left(\frac{1}{\lambda }\right)$$

where $D_{KL}[P\left(y\right|\lambda )\parallel P_{AC}\left(y\right|x)]$ is the Kullback–Leibler divergence from $P\left(y\right|\lambda )$ to $P_{AC}\left(y\right|x)$,

$$D_{KL}[P\left(y\right|\lambda )\parallel P_{AC}\left(y\right|x)]=\sum _{y=0}^{\infty }P\left(y\right|\lambda )log\frac{P\left(y\right|\lambda )}{P_{AC}\left(y\right|x)}$$

The expected divergence (in bits) can be well-approximated (up to order $O\left({\lambda }^{-3}\right)$) by [11]:

$$\mathcal{E}\left(\lambda \right)\approx \frac{1}{2}-\frac{1}{12\lambda }\left(1-\frac{1}{2}\right)-\frac{1}{24{\lambda }^2}\left(1-\frac{1}{2^2}\right)$$

(6)

2.3. $P_{AC}\left(y\right|x)$ vs. Maximum Likelihood

In this section we will briefly recall information theoretic analysis of the maximum likelihood estimate $P_{ML}\left(y\right|x)$ in place of $P_{AC}\left(y\right|x)$ [12]. First note that Poisson distribution $P\left(y\right|\lambda )$ is only defined for positive λ. In the case of observing zero count $x=0$, we cannot directly use the “maximum likelihood estimate” $P\left(y\right|0)$. One option for dealing with zero observed counts is to allow for some form of model regularization, e.g., infer a Poisson model $P\left(y\right|ϵ)$, for some small $ϵ>0$. In other words, if a count $x\ge 1$ is observed, follow the standard maximum likelihood procedure and infer $P_{ML}\left(y\right|x)=P\left(y\right|x)$ as the Poisson model; if a zero count is observed, $x=0$, infer $P_{ML}\left(y\right|0)=P\left(y\right|ϵ)$ for some fixed $ϵ\in (0,1]$. This is the route taken in [12] and adopted in this paper. Only a minimum amount of necessary regularization due to zero observed counts is employed in the otherwise straightforward ML approach.

Theorem 2 [12] Consider an underlying Poisson distribution $P(·|\lambda )$ parameterized by some $\lambda >0$ and a regularization constant $ϵ\in (0,1]$. The expected divergence in bits $\mathsf{{\rm Y}}(\lambda ,ϵ)$ between the true Poisson source and its (regularized) maximum likelihood estimate based on a single observation,

$$\mathsf{{\rm Y}}(\lambda ,ϵ)=E_{P\left(x\right|\lambda )}\left[D_{KL}[P\left(y\right|\lambda )\parallel P_{ML}\left(y\right|x)]\right]$$

is equal to

$$\mathsf{{\rm Y}}(\lambda ,ϵ)=\lambda \left({log}_2\lambda -\sum _{x=1}^{\infty }P\left(x\right|\lambda ){log}_{2}x\right)+e^{-\lambda }(ϵ-\lambda {log}_2ϵ)$$

(7)

Note that the expected divergence $\mathsf{{\rm Y}}(\lambda ,ϵ)$ can get prohibitively large when regularizing with small $ϵ>0$. As an illustration, in Figure 1 we show expected divergence $\mathsf{{\rm Y}}(\lambda ,ϵ=1)$ of the ML estimation (zero count regularized with $ϵ=1$) for a range of mean parameter values λ of the underlying Poisson source (solid line). Also shown is the expected divergence $\mathcal{E}\left(\lambda \right)$ of $P_{AC}\left(y\right|x)$ (dashed line). Except for very small Poisson source rates λ, $P_{AC}\left(y\right|x)$ is clearly benefitting from the stabilizing effect of Bayesian averaging, given the extremely small sample size.

Figure 1. Expected divergence (in bits) $\mathsf{{\rm Y}}(\lambda ,ϵ=1)$ of the ML estimation (zero count regularized with $ϵ=1$) (solid line). Also shown is the expected divergence $\mathcal{E}\left(\lambda \right)$ of $P_{AC}\left(y\right|x)$ (dashed line).

Figure 1. Expected divergence (in bits) $\mathsf{{\rm Y}}(\lambda ,ϵ=1)$ of the ML estimation (zero count regularized with $ϵ=1$) (solid line). Also shown is the expected divergence $\mathcal{E}\left(\lambda \right)$ of $P_{AC}\left(y\right|x)$ (dashed line).

3. Generalized $P_{AC}\left(y\right|x)$ with Gamma Prior

In this section we will generalize $P_{AC}\left(y\right|x)$ through the use of (conjugate) gamma prior

$$P\left(\lambda \right|\alpha ,\beta )=\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{\lambda }^{\alpha -1}{e}^{-\beta \lambda }$$

on the Poisson mean parameter λ. The positive parameters $\alpha ,\beta$ determine the overall shape of the prior. Given a single observation x, the posterior

$$P\left(\lambda \right|x,\alpha ,\beta )=\frac{P\left(x\right|\lambda \left)P\right(\lambda |\alpha ,\beta )}{{\int }_0^{\infty }P\left(x\right|\lambda )P\left(\lambda \right|\alpha ,\beta )\text{d}\lambda }$$

is the gamma distribution with parameters $\alpha +x$ and $\beta +1$,

$$P\left(\lambda \right|x,\alpha ,\beta )=\frac{{(\beta +1)}^{\alpha +x}}{\Gamma (\alpha +x)}{\lambda }^{\alpha +x-1}{e}^{-(\beta +1)\lambda }$$

The mean of $P\left(\lambda \right|x,\alpha ,\beta )$ is equal to $(\alpha +x)/(\beta +1)$. A loose intuitive interpretation of the prior parameters $\alpha ,\beta$ (assuming they are integers) is that prior to seeing the current data (in our case only one observation (count) x), we have seen β “observations”, $x_1^{\prime },x_2^{\prime },...,x_{\beta }^{\prime }$, with the total cumulative count $\alpha =x_1^{\prime }+x_2^{\prime }+...+x_{\beta }^{\prime }$. Hence the mean parameter estimate would shift from x (ML estimation corresponding to $\alpha ,\beta \to 0$) to $(x_1^{\prime }+x_2^{\prime }+...+x_{\beta }^{\prime }+x)/(\beta +1)$.

As in the case of $P_{AC}\left(y\right|x)$, having observed a count x, we build a predictive distribution over future counts y by integrating out the mean parameter λ with respect to the posterior $P\left(\lambda \right|x,\alpha ,\beta )$,

$$\begin{array}{ccc}\hfill P_G\left(y\right|x,\alpha ,\beta )& =& {\int }_0^{\infty }P\left(y\right|\lambda )P\left(\lambda \right|x,\alpha ,\beta )\text{d}\lambda \hfill \\ & =& \frac{{(\beta +1)}^{\alpha +x}}{\Gamma (\alpha +x)}\frac{1}{y!}{\int }_0^{\infty }{\lambda }^{\alpha +x+y-1}{e}^{-(\beta +2)\lambda }\text{d}\lambda \hfill \end{array}$$

(8)

From normalization of the gamma distribution we get

$${\int }_0^{\infty }{\lambda }^{a-1}{e}^{-b\lambda }\text{d}\lambda =\frac{\Gamma \left(a\right)}{b^a}$$

and so

$${\int }_0^{\infty }{\lambda }^{\alpha +x+y-1}{e}^{-(\beta +2)\lambda }\text{d}\lambda =\frac{\Gamma (x+y+\alpha )}{{(\beta +2)}^{x+y+\alpha }}$$

leading to

$$P_G\left(y\right|x,\alpha ,\beta )=\frac{1}{y!}\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}$$

(9)

It can be easily verified that the original $P_{AC}\left(y\right|x)$ is obtained as a special case of $P_G\left(y\right|x,\alpha ,\beta )$ when $\alpha =1$ and $\beta \to 0$. If Jeffrey’s prior were used instead of the flat prior in $P_{AC}\left(y\right|x)$, we would obtain $P_G\left(y\right|x,\alpha ,\beta )$ with $\alpha =1/2$ and $\beta \to 0$ etc.

If α is an integer, we have

$$P_G\left(y\right|x,\alpha ,\beta )={\left(\frac{1+\beta }{2+\beta }\right)}^{x^{\prime }+1}{\left(\frac{1}{2+\beta }\right)}^y\left(\genfrac{}{}{0pt}{}{x^{\prime }+y}{y}\right)$$

(10)

where $x^{\prime }=x+\alpha -1$ is the observed count including prior observations. This expression generalizes $P_{AC}\left(y\right|x)$ (5),

$$P_{AC}\left(y\right|x)={\left(\frac{1}{2}\right)}^{x+1}{\left(\frac{1}{2}\right)}^y\left(\genfrac{}{}{0pt}{}{x+y}{y}\right)$$

While $P_G\left(y\right|x,\alpha ,\beta )$ (9) can be used with any appropriate setting of $\alpha ,\beta$ (e.g., given a prior knowledge of the range of counts one may reasonably expect), in this contribution we concentrate on using the gamma prior to mitigate for the unrealistic equal weighting of all $\lambda >0$ in the flat prior behind $P_{AC}\left(y\right|x)$. Indeed, the observed counts are typically bounded by the nature of the problem and one can represent this through setting $\alpha =1$ and varying $\beta >0$ in the gamma prior $P\left(\lambda \right|\alpha ,\beta )$ underlying $P_G\left(y\right|x,\alpha ,\beta )$. Some examples of such priors are shown in Figure 2. Decreasing β leads to weaker emphasis on low λ, eventually recovering the flat (improper) prior for $\beta =0$.

Figure 2. Gamma prior $P\left(\lambda \right|\alpha =1,\beta )$. Shown are the priors for three possible values of parameter β, $\beta \in \{1,0.1,0.05\}$.

Figure 2. Gamma prior $P\left(\lambda \right|\alpha =1,\beta )$. Shown are the priors for three possible values of parameter β, $\beta \in \{1,0.1,0.05\}$.

In Section 2.3 maximum likelihood estimation was regularized at zero count by imposing a non-zero “count” ϵ instead of the observed zero one. The generalized form of $P_{AC}\left(y\right|x)$, $P_R\left(y\right|x,\beta )=P_G\left(y\right|x,\alpha =1,\beta )$ can be also viewed as an alternative “soft” form of regularization of the maximum likelihood approach at zero counts.

Parameter β in the Gamma prior

$$P\left(\lambda \right|\alpha =1,\beta )=\beta e^{-\beta \lambda }$$

can be set in a data driven manner, e.g., using the following strategy: Given the observed count x, we require that the area up to $x+1$ covered by the prior is equal to θ, for some threshold $\theta \in (0,1)$ (e.g., $\theta =1/4$). In other words, $F(x+1|\beta )=\theta$, where $F\left(\lambda \right|\beta )=1-e^{-\beta \lambda }$ is the cumulative distribution function of $P\left(\lambda \right|\alpha =1,\beta )$. This leads to

$$\beta \left(x\right)=-\frac{ln(1-\theta )}{x+1}$$

(11)

For zero observed count $x=0$, $\beta \left(0\right)=-ln(1-\theta )$ and the prior gets more concentrated on smaller values of λ as likely candidates for the mean count of the underlying Poisson source. With increasing count values $x>0$ the parameter $\beta \left(x\right)$ decreases to 0 and the prior gradually approaches the flat prior of $P_{AC}\left(y\right|x)$.

Finally, we contrast $P_G\left(y\right|x,\alpha ,\beta )$ with the negative binomial distribution

$$P_{NB}\left(y\right|r,q)=\frac{1}{y!}\frac{\Gamma (r+y)}{\Gamma \left(r\right)}{q}^r{(1-q)}^y$$

(12)

with parameters $r>0$ and $q\in [0,1]$. One interpretation of the negative binomial distribution $P_{NB}\left(y\right|r,q)$ is that it corresponds to a Gamma–Poisson mixture that one obtains by imposing a Gamma prior $P\left(\lambda \right|r,(1-q)/q)$ on the mean count parameter λ of the Poisson distribution $P\left(y\right|\lambda )$ and integrating out λ. In our context it is natural to identify r and $(1-q)/q$ with hyperparameters α and β used in $P_G\left(y\right|x,\alpha ,\beta )$. It follows that $q={(\beta +1)}^{-1}$. Hence, we rewrite Equation (12) as

$$P_{NB}\left(y\right|\alpha ,{(\beta +1)}^{-1})=\frac{1}{y!}\frac{\Gamma (\alpha +y)}{\Gamma \left(\alpha \right)}\frac{{\beta }^y}{{(\beta +1)}^{\alpha +y}}$$

(13)

Direct comparison of (13) with (9) leads to an intuitive insight: The β prior measurements of total count α introduced by the gamma prior $P\left(\lambda \right|\alpha ,\beta )$ are in the case of $P_G\left(y\right|x,\alpha ,\beta )$ extended with a single observation x, resulting in $\beta +1$ observations of total count $\alpha +x$. This can be represented by

$$P_{NB}\left(y\right|\alpha +x,{(\beta +2)}^{-1})=\frac{1}{y!}\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^y}{{(\beta +2)}^{x+y+\alpha }}$$

(14)

It follows that

$$\frac{P_G\left(y\right|x,\alpha ,\beta )}{P_{NB}\left(y\right|\alpha +x,{(\beta +2)}^{-1})}={(\beta +1)}^{x+\alpha -y}.$$

Bayesian averaging in $P_G\left(y\right|x,\alpha ,\beta )$ with respect to the posterior over λ, given a count x, differs from the corresponding negative binomial distribution $P_{NB}\left(y\right|\alpha +x,{(\beta +2)}^{-1})$ by the factor ${(\beta +1)}^{x+\alpha -y}$ that depends on the difference between the prior+observed count $\alpha +x$ and y.

4. First and Second Moments of the Generalized $P_{AC}\left(y\right|x)$

In [11] we showed that $P_{AC}\left(y\right|x)$ and the underlying Poisson distribution are quite similar in their nature: for any (integer) mean rate $\lambda \ge 1$, the Poisson distribution $P(·|\lambda )$ has two neighboring modes located at λ and $\lambda -1$, with $P\left(\lambda \right|\lambda )=P(\lambda -1\left|\lambda \right)$. Analogously, given a count $x\ge 1$, $P_{AC}(·|x)$ has two neighboring modes, one located at x, the other at $x-1$, with $P_{AC}\left(x\right|x)=P_{AC}(x-1|x)$. As in Poisson distribution, the values of $P_{AC}\left(y\right|x)$ decrease as one moves away from the modes in both directions. In this section we derive the first two moments of the generalized $P_{AC}\left(y\right|x)$, $P_G\left(y\right|x,\alpha ,\beta )$. As a special case, we will show that as a result of Bayesian averaging, the variance of $P_{AC}\left(y\right|x)$ is double that of the underlying (unobserved) Poisson distribution.

Theorem 3 Consider a non-negative integer x and the associated generalized model $P_G\left(y\right|x,\alpha ,\beta )$. Then,

$$E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]=\frac{x+\alpha }{\beta +1},Var\left[y\right]=\frac{\beta +2}{\beta +1}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]$$

Proof: Let us evaluate

$$\begin{array}{ccc}\hfill E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]& =& \sum _{y=0}^{\infty }\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}\frac{1}{y!}y\hfill \\ & =& \sum _{y=1}^{\infty }\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}\frac{1}{(y-1)!}\hfill \\ & =& \sum _{y^{\prime }=0}^{\infty }\frac{\Gamma (x+y^{\prime }+1+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y^{\prime }+1+\alpha }}\frac{1}{y^{\prime }!}\hfill \\ & =& \sum _{y^{\prime }=0}^{\infty }\frac{\Gamma (x+y^{\prime }+\alpha )·(x+y^{\prime }+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{(\beta +2)·{(\beta +2)}^{x+y^{\prime }+\alpha }}\frac{1}{y^{\prime }!}\hfill \end{array}$$

(15)

In the third equality we have used substitution $y^{\prime }=y-1$ and the last equality follows from $\Gamma (z+1)=z·\Gamma \left(z\right)$. By (15),

$$\begin{array}{ccc}\hfill E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]& =& \sum _{y=0}^{\infty }{P}_G\left(y\right|x,\alpha ,\beta )\frac{x+\alpha +y}{\beta +2}\hfill \end{array}$$

(16)

$$\begin{array}{ccc}\phantom{\hfill E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]}& =& \frac{x+\alpha }{\beta +2}+\frac{1}{\beta +2}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]\hfill \end{array}$$

(17)

Solving (17) we obtain

$$E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]=\frac{x+\alpha }{\beta +1}$$

(18)

For the variance of $P_G\left(y\right|x,\alpha ,\beta )$ we have

$$\begin{array}{ccc}\hfill Va{r}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]& =& E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]-{\left(E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]\right)}^2\hfill \end{array}$$

(19)

Now,

$$\begin{array}{ccc}\hfill E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]& =& \sum _{y=0}^{\infty }\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}\frac{1}{y!}{y}^2\hfill \\ & =& \sum _{y=1}^{\infty }\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}\frac{1}{y!}{y}^2\hfill \\ & =& \sum _{y=1}^{\infty }\frac{\Gamma (x+y+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y+\alpha }}\frac{1}{(y-1)!}y\hfill \\ & =& \sum _{y^{\prime }=0}^{\infty }\frac{\Gamma (x+y^{\prime }+1+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{{(\beta +2)}^{x+y^{\prime }+1+\alpha }}\frac{1}{y^{\prime }!}(y^{\prime }+1)\hfill \\ & =& \sum _{y^{\prime }=0}^{\infty }\frac{(x+y^{\prime }+\alpha )\Gamma (x+y^{\prime }+\alpha )}{\Gamma (x+\alpha )}\frac{{(\beta +1)}^{x+\alpha }}{(\beta +2){(\beta +2)}^{x+y^{\prime }+\alpha }}\frac{1}{y^{\prime }!}(y^{\prime }+1)\hfill \\ & =& \frac{1}{\beta +2}\sum _{y^{\prime }=0}^{\infty }{P}_G\left(y^{\prime }\right|x,\alpha ,\beta )\left[(x+y^{\prime }+\alpha )(y^{\prime }+1)\right]\hfill \\ & =& \frac{1}{\beta +2}\sum _{y^{\prime }=0}^{\infty }{P}_G\left(y^{\prime }\right|x,\alpha ,\beta )[x+y^{\prime }+\alpha ]\hfill \\ & +& \frac{1}{\beta +2}\sum _{y^{\prime }=0}^{\infty }{P}_G\left(y^{\prime }\right|x,\alpha ,\beta )[y^{\prime }(x+\alpha )+y^{\prime 2}]\hfill \end{array}$$

(20)

Using (16), (18) and (20), we obtain

$$\begin{array}{ccc}\hfill E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]& =& E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]+\frac{x+\alpha }{\beta +2}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]+\frac{1}{\beta +2}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]\hfill \\ \hfill \\ & =& \frac{x+\alpha }{\beta +1}\left(1+\frac{x+\alpha }{\beta +2}\right)+\frac{1}{\beta +2}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]\hfill \end{array}$$

(21)

which can be solved as

$$E_{P_G\left(y\right|x,\alpha ,\beta )}\left[y^2\right]=\frac{(x+\alpha )(x+\alpha +\beta +2)}{{(\beta +1)}^2}$$

(22)

Plugging (22) into (19) we obtain

$$Va{r}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]=\frac{(x+\alpha )(\beta +2)}{{(\beta +1)}^2}=\frac{\beta +2}{\beta +1}{E}_{P_G\left(y\right|x,\alpha ,\beta )}\left[y\right]$$

                              ☐

Given an observation x, the maximum likelihood estimate of the underlying Poisson distribution is the Poisson distribution with mean x,

$$P\left(y\right|x)=e^{-x}\frac{x^y}{y!}$$

After observing x, the mean of the maximum likelihood and $P_{AC}(·|x)$ estimates is x and $x+1$, respectively. Hence, Bayesian averaging in $P_{AC}(·|x)$ induced by the flat improper prior over the mean rate λ results in increased expected value $x+1$ of the next count from the same underlying source, given that the current count x. However, a much more marked consequence of using the flat prior can be seen in the variance of $P_{AC}(·|x)$: while variance of the maximum likelihood is x, it is $2(x+1)$ in $P_{AC}(·|x)$.

Theorem 3 illustrates the role of more concentrated prior over λ on the generalized model. The mean expected count, after seeing x, is equal to the mean of the posterior $P\left(\lambda \right|x,\alpha ,\beta )$ over λ, namely $(\alpha +x)/(\beta +1)$. As explained earlier, observed single count x with prior β counts of cumulative value α results in $\beta +1$ counts of cumulative value $\alpha +x$. Hence the mean count per observation is $(\alpha +x)/(\beta +1)$. As with Poisson distribution, the variance of the generalized model is closely related to its mean and approaches the mean with increasing number of prior counts β.

As for the soft regularization $P_R\left(y\right|x,\beta )=P_G\left(y\right|x,\alpha =1,\beta )$, its mean is, as expected, biased towards values smaller than the observed count x, provided $\beta >1/x$. Increased values of β result in smaller variance of $P_R\left(y\right|x,\beta )$. But how do such prior parameter modifications manifest themselves in terms of accuracy of estimation of the underlying source? This question is investigated in the next section.

5. Expected Divergence of the Generalized $P_{AC}\left(y\right|x)$ from the True Underlying Poisson Distribution

Consider an underlying Poisson source $P\left(x\right|\lambda )$ generating counts x. In this section we would like to quantify the average divergence

$${\mathcal{E}}_G(\lambda ;\beta )=E_{P\left(x\right|\lambda )}\left[D_{KL}[P\left(y\right|\lambda )\parallel P_R\left(y\right|x,\beta )]\right]$$

(23)

of the corresponding generalized $P_{AC}\left(y\right|x)$, $P_R\left(y\right|x,\beta )=P_G\left(y\right|x,\alpha =1,\beta )$ (“softly” regularized ML), from the truth $P\left(y\right|\lambda )$, if we repeatedly generated a “representative” count x from $P\left(x\right|\lambda )$. The same question was considered in the context of maximum likelihood estimation in Section 2.3. In particular, we are interested in specifying under what circumstances is the generalized form of $P_{AC}\left(y\right|x)$, $P_R\left(y\right|x,\beta )=P_G\left(y\right|x,\alpha =1,\beta )$, preferable to the original $P_{AC}\left(y\right|x)=P_G\left(y\right|x,\alpha =1,\beta \to 0)$ and how it fares with the maximum likelihood estimation $P_{ML}\left(y\right|x)$ of Section 2.3.

Theorem 4 Consider an underlying Poisson distribution $P(·|\lambda )$ parameterized by some $\lambda >0$. Then for $\beta \ge 0$,

$${\mathcal{E}}_G(\lambda ;\beta )={log}_2\left(\frac{\beta +2}{\beta +1}\right)-\frac{1}{2}+\lambda \left[2{log}_2\left(\frac{\beta +2}{2}\right)-{log}_2(\beta +1)\right]+O\left({\lambda }^{-1}\right)$$

(24)

A higher order approximation (up to order ${\lambda }^3$) reads:

$$\begin{array}{ccc}\hfill {\mathcal{E}}_G(\lambda ;\beta )& =& {log}_2\left(\frac{\beta +2}{\beta +1}\right)-\frac{1}{2}+\lambda \left[2{log}_2\left(\frac{\beta +2}{2}\right)-{log}_2(\beta +1)\right]\hfill \\ & -& \frac{1}{12\lambda }\left(1-\frac{1}{2}\right)-\frac{1}{24{\lambda }^2}\left(1-\frac{1}{2^2}\right)-\frac{19}{360{\lambda }^3}\left(1-\frac{1}{2^3}\right)+O\left({\lambda }^{-4}\right)\hfill \end{array}$$

(25)

Proof: Let us first express the divergence $D_{\beta }(\lambda ,x)=D_{KL}[P\left(y\right|\lambda )\parallel P_R\left(y\right|x,\beta )]$. We have

$$D_{\beta }(\lambda ,x)=-H\left[P\left(y\right|\lambda )\right]-E_{P\left(y\right|\lambda )}[log{P}_R\left(y\right|x,\beta )]$$

where $H\left[P\left(y\right|\lambda )\right]=-E_{P\left(y\right|\lambda )}[logP\left(y\right|\lambda )]$ is the entropy of the source $P\left(y\right|\lambda )$ and

$$\begin{array}{ccc}\hfill E_{P\left(y\right|\lambda )}[log{P}_R\left(y\right|x,\beta )]& =& -logx!\hfill \\ & & -E_{P\left(y\right|\lambda )}\left[y\right]log(\beta +2)-(x+1)log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & -E_{P\left(y\right|\lambda )}[logy!]+E_{P\left(y\right|\lambda )}[log(x+y)!]\hfill \end{array}$$

Denoting (for integer $d\ge 0$) $E_{P\left(y\right|\lambda )}[log(y+d)!]$ by $F(\lambda ,d)$, we write

$$\begin{array}{ccc}\hfill D_{\beta }(\lambda ,x)& =& -H\left[P\right(y\left|\lambda \right)]+logx!\hfill \\ & & +\lambda log(\beta +2)+(x+1)log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & +F(\lambda ,0)-F(\lambda ,x)\hfill \end{array}$$

We are now ready to calculate the expectation ${\mathcal{E}}_G(\lambda ;\beta )=E_{P\left(x\right|\lambda )}\left[D_{\beta }(\lambda ,x)\right]$.

$$\begin{array}{ccc}\hfill {\mathcal{E}}_G(\lambda ;\beta )& =& -H\left[P\right(y\left|\lambda \right)]+F(\lambda ,0)\hfill \\ & & +\lambda log(\beta +2)+(\lambda +1)log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & +F(\lambda ,0)-E_{P\left(x\right|\lambda )}\left[F(\lambda ,x)\right]\hfill \end{array}$$

We have proved in [11] that $E_{P\left(x\right|\lambda )}\left[F(\lambda ,x)\right]=F(2\lambda ,0)$, and so

$$\begin{array}{ccc}\hfill {\mathcal{E}}_G(\lambda ;\beta )& =& -H\left[P\left(y\right|\lambda )\right]+log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & +\lambda log\left(\frac{{(\beta +2)}^2}{\beta +1}\right)\hfill \\ & & +2F(\lambda ,0)-F(2\lambda ,0)\hfill \end{array}$$

Since

$$\begin{array}{ccc}\hfill -H\left[P\right(y\left|\lambda \right)]& =& E_{P\left(y\right|\lambda )}[logP\left(y\right|\lambda )]\hfill \\ & =& -\lambda loge+E_{P\left(y\right|\lambda )}\left[y\right]log\lambda -E_{P\left(y\right|\lambda )}[logy!]\hfill \\ & =& -\lambda loge+\lambda log\lambda -F(\lambda ,0)\hfill \end{array}$$

(26)

we have

$$\begin{array}{ccc}\hfill {\mathcal{E}}_G(\lambda ;\beta )& =& log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & +\lambda \left[log\lambda +log\left(\frac{{(\beta +2)}^2}{\beta +1}\right)-loge\right]\hfill \\ & & +F(\lambda ,0)-F(2\lambda ,0)\hfill \end{array}$$

(27)

Using entropy approximation (see [11]), one obtains

$$F(\lambda ,0)=\lambda (log\lambda -loge)+\frac{1}{2}log\left(2\pi e\lambda \right)+O\left({\lambda }^{-1}\right)$$

leading to (in log base 2)

$$F(\lambda ,0)-F(2\lambda ,0)=-\frac{1}{2}+\lambda ({log}_{2}e-{log}_2\lambda -2)+O\left({\lambda }^{-1}\right)$$

Finally,

$$\begin{array}{ccc}\hfill {\mathcal{E}}_G(\lambda ;\beta )& =& {log}_2\left(\frac{\beta +2}{\beta +1}\right)-\frac{1}{2}\hfill \\ & & +\lambda \left[{log}_2\left(\frac{{(\beta +2)}^2}{\beta +1}\right)-2\right]+O\left({\lambda }^{-1}\right)\hfill \end{array}$$

which is equivalent to (24).

The higher order expression (25) is simply obtained by using higher order approximation to $F(\lambda ,0)-F(2\lambda ,0)$.

☐

Note that for $\beta \to 0$ we recover our original result [11] that the expected divergence $\mathcal{E}\left(\lambda \right)$ of the original $P_{AC}\left(y\right|x)$ from the “truth” $P\left(y\right|\lambda )$ is (up to terms of order ${\lambda }^{-1})$ never greater than 1/2 bit. The soft regularization in $P_R\left(y\right|x,\beta )$ (using prior $P\left(\lambda \right|\alpha =1,\beta )$ with $\beta >0$) can result in larger expected divergence from the underlying source than is the case for $P_{AC}\left(y\right|x)$ (using improper flat prior over λ). Moreover, (unlike in $P_{AC}\left(y\right|x)$) such a regularization causes linear divergence of ${\mathcal{E}}_G(\lambda ;\beta )$ for large λ. The next theorem specifies for which underlying Poisson sources the soft regularization approach of $P_R\left(y\right|x,\beta )$ is preferable to the original $P_{AC}\left(y\right|x)$.

Theorem 5 For Poisson sources with mean rates

$$\lambda <\kappa \left(\beta \right)=\frac{log\left(1+\frac{\beta }{\beta +2}\right)}{log\left(1+\frac{{\beta }^2}{4(\beta +1)}\right)}$$

(28)

it holds $\mathcal{E}\left(\lambda \right)>{\mathcal{E}}_G(\lambda ;\beta )$ and hence $P_R\left(y\right|x,\beta )$ is on average guaranteed to approximate (in the Kullback–Leibler divergence sense) the underlying source better than the original $P_{AC}\left(y\right|x)$.

Proof: It was shown in [11] that for the original $P_{AC}\left(y\right|x)$,

$$\begin{array}{ccc}\hfill \mathcal{E}\left(\lambda \right)& =& \lambda (log\lambda -loge+2log2)+log2\hfill \\ & & +F(\lambda ,0)-F(2\lambda ,0).\hfill \end{array}$$

(29)

From (27) and (29) we have that the difference between the expected divergences of the original and generalized forms of $P_{AC}\left(y\right|x)$ is

$$\begin{array}{ccc}\hfill \mathcal{E}\left(\lambda \right)-{\mathcal{E}}_G(\lambda ;\beta )& =& log2-log\left(\frac{\beta +2}{\beta +1}\right)\hfill \\ & & +\lambda \left[2log2-log\left(\frac{{(\beta +2)}^2}{\beta +1}\right)\right]\hfill \\ & =& log\frac{2(\beta +1)}{\beta +2}\hfill \\ & & +\lambda log\frac{4(\beta +1)}{{(\beta +2)}^2}\hfill \end{array}$$

(30)

The result follows from solving for $\mathcal{E}\left(\lambda \right)>{\mathcal{E}}_G(\lambda ;\beta )$.

The graph (in log-log scale) of $\kappa \left(\beta \right)$ is shown in Figure 3. An alternative way of data-driven setting of parameter β is suggested by the fact that $\kappa \left(\beta \right)$ is lower bounded by ${\beta }^{-1}$. If the experimental setting is such that most counts are expected not to exceed some $x_{max}$, β can be set to $\beta =1/x_{max}$, so that $P_R\left(y\right|x,\beta )$ is preferable to $P_{AC}\left(y\right|x)$.

In Figure 4 we present the expected divergences ${\mathcal{E}}_G(\lambda ;\beta )$ (solid line) and $\mathcal{E}\left(\lambda \right)$ (dashed line) for $\beta =0.2$ (left) and $\beta =0.01$ (right). As expected, for underlying sources with small mean counts λ the advantage of using the regularized form $P_R\left(y\right|x,\beta )$ (as opposed to the original $P_{AC}\left(y\right|x)$) is more pronounced. However, for larger λ there is a heavy price to be paid in terms of inaccurate modelling by $P_R\left(y\right|x,\beta )$.

Figure 3. Graph of $\kappa \left(\beta \right)$. For Poisson sources with mean rates $\lambda <\kappa \left(\beta \right)$, $\mathcal{E}\left(\lambda \right)>{\mathcal{E}}_G(\lambda ;\beta )$ and hence $P_R\left(y\right|x,\beta )$ is on average guaranteed to approximate the underlying source better than the original $P_{AC}\left(y\right|x)$.

Figure 3. Graph of $\kappa \left(\beta \right)$. For Poisson sources with mean rates $\lambda <\kappa \left(\beta \right)$, $\mathcal{E}\left(\lambda \right)>{\mathcal{E}}_G(\lambda ;\beta )$ and hence $P_R\left(y\right|x,\beta )$ is on average guaranteed to approximate the underlying source better than the original $P_{AC}\left(y\right|x)$.

Figure 4. Expected divergences ${\mathcal{E}}_G(\lambda ;\beta )$ (solid line) and $\mathcal{E}\left(\lambda \right)$ (dashed line) for $\beta =0.2$ (left) and $\beta =0.01$ (right).

Figure 4. Expected divergences ${\mathcal{E}}_G(\lambda ;\beta )$ (solid line) and $\mathcal{E}\left(\lambda \right)$ (dashed line) for $\beta =0.2$ (left) and $\beta =0.01$ (right).

6. Empirical Investigations

To investigate potential value of the more sophisticated Bayesian approach in the original and the generalized Audic–Claverie frameworks (Section 2.1 and Section 3, respectively) against the baseline of simple (regularized) maximum likelihood estimation (Section 2.3), we conducted a series of simple illustrative experiments. In the generalized Audic–Claverie framework developed in this study, we used the two schemes for setting the regularization parameter β suggested in Section 3 and Section 5. In the regularized maximum likelihood approach $P_{ML}\left(y\right|x)$ we set $ϵ=1$. From Figure 1, it appears that the biggest difference between the expected divergences from the true underlying Poisson source $P\left(x\right|\lambda )$ to the original $P_{AC}(·|x)$ and the maximum likelihood estimate occurs for small mean rates λ roughly around $\lambda =5$. We therefore run the experiments with $\lambda =5$.

For illustration purposes, we follow the data generation mechanism used in [13] to compare methods for distinguishing between differential expression of genes associated with two treatment regimes. We stress that in no way we suggest that our experiments have strong relevance for bioinformatics, nor do we claim that the framework of [13] is the best test bed for assessing differential gene expression detection algorithms. We use the framework of [13] merely to illustrate whether the sophistication of the Bayesian approach (as opposed to simple (regularized) maximum likelihood) can bring benefits in a practical situation with low-count data.

Gene counts are simulated across the two treatment groups $T_1$ and $T_2$. The tests are assessed by comparing false positive and true positive rates. In each experiment 10,000 gene pair counts $(x_{1,j},x_{2,j})$, $j=1,2,...,10,000$, were produced, counts $x_{1,j}$ and $x_{2,j}$ associated with regimes $T_1$ and $T_2$, respectively. As specified above, the sampling rate for $T_1$ was fixed at ${\lambda }_1=5$ throughout the experiment. We varied the mean ${log}_2$ fold change (LFC) between $T_1$ and $T_2$ from −2 to 2. Each gene pair count $(x_{1,j},x_{2,j})$, $j=1,2,...,10,000$, was obtained through a generative process specified in [13] and described in detail in Appendix A.

Having generated the gene pair counts, we used methods considered in this study to make a decision for each $j=1,2,...,10,000$, whether the counts $x_{1,j},x_{2,j}$ originated from the same underlying source, i.e., whether when generating $y_{1,j}$ and $y_{2,j}$, the mean rates in the two regimes $T_1$ and $T_2$ were identical ($LF{C}_j=0$). Given the “test distribution” $Q\left(y\right|x)$ and a confidence level $\vartheta \in [0,1]$, we guess that $x_{1,j},x_{2,j}$ originated from the same source if the $(1-\vartheta )$-quantile around the mean of $Q\left(y\right|x_{1,j})$ contains $x_{2,j}$ and vice-versa, i.e., if the $(1-\vartheta )$-quantile around the mean of $Q\left(y\right|x_{2,j})$ contains $x_{1,j}$. In place of $Q\left(y\right|x)$ we used $P_{AC}\left(y\right|x)$, its regularized form $P_R\left(y\right|x,\beta )$ and the regularized maximum likelihood estimate $P_{ML}\left(y\right|x)$ with $ϵ=1$.

For a given confidence level $\vartheta \in [0,1]$ and test statistic $Q\left(y\right|x)$ we calculate the false positive rate (type I error rate) as the proportion of times a gene count pair $(x_{1,j},x_{2,j})$ was declared to have originated from two different underlying sources (differentially expressed gene) when in fact $LF{C}_j$ was zero. The true positive rate (statistical power) was determined as the proportion of times a gene was correctly declared differentially expressed - $(x_{1,j},x_{2,j})$ declared to have originated come from two different underlying sources and $LF{C}_j\ne 0$.

Plot of false positive rate vs. true positive rate obtained for different values of ϑ constitutes a receiver operating characteristic (ROC) curve. If the ROC curve for one test distribution is always above another, this suggests its superiority in classifying genes as differentially expressed. Trivial classification of genes as differentially expressed using a completely random guess would yield the identity (diagonal) ROC curve. ROC curves for the maximum likelihood method ($ϵ=1$, red dashed line) and the soft regularization model $P_R\left(y\right|x,\beta )$, $\beta =1/50,1/100$ (solid lines) are plotted in Figure 5. Not surprisingly, the Bayesian approach (solid lines) outperforms the penalized maximum likelihood one (red dashed line). However, the original $P_{AC}\left(y\right|x)$ ($\beta =0$, black line) and the soft regularization model (color solid lines) achieve almost identical performances. In this challenging setting (single observations at low mean rate with additional noise), the scheme for setting the regularization parameter β suggested in Section 5 has little effect on the resulting classification performance. We also ran experiments to test the “dynamic” scheme for setting β introduced in Section 3, but no significant performance improvements were achieved.

Figure 5. ROC curves for test distributions $P_{AC}\left(y\right|x)=P_R\left(y\right|x,\beta \to 0)$ (solid black line), $P_R\left(y\right|x,\beta =1/100)$ (solid blue line), $P_R\left(y\right|x,\beta =1/50)$ (solid green line) and $P_{ML}\left(y\right|x)$ with $ϵ=1$ (dashed red line). Mean rate of the underlying Poisson source was fixed at $\lambda =5$.

Figure 5. ROC curves for test distributions $P_{AC}\left(y\right|x)=P_R\left(y\right|x,\beta \to 0)$ (solid black line), $P_R\left(y\right|x,\beta =1/100)$ (solid blue line), $P_R\left(y\right|x,\beta =1/50)$ (solid green line) and $P_{ML}\left(y\right|x)$ with $ϵ=1$ (dashed red line). Mean rate of the underlying Poisson source was fixed at $\lambda =5$.

Finally, we devised yet another scheme for determining the hyper-parameters α and β of the prior $P\left(\lambda \right|\alpha ,\beta )$ from the data. In the spirit of type II maximum likelihood, we find the most likely values of $\alpha ,\beta$, given the observed counts $\mathcal{C}=\{x_1,x_2,...,x_n\}$, using $P\left(\mathcal{C}\right|\alpha ,\beta )={\prod }_{i=1}^{n}P\left(x_i\right|\alpha ,\beta ),$ where

$$P\left(x_i\right|\alpha ,\beta )={\int }_0^{\infty }P\left(x_i\right|\lambda )p\left(\lambda \right|\alpha ,\beta )\text{d}\lambda$$

(31)

Using this method, we first optimize the prior hyperparameters on the observed data. The “optimized” prior $P\left(x_i\right|{\alpha }_{*},{\beta }_{*})$ now reflects the possible ranges of mean counts λ one can expect given the data. We then repeated the experiments using the generalized model $P_G\left(y\right|x,{\alpha }_{*},{\beta }_{*})$ derived from the optimized prior. In this way we can assess to what degree the relatively minor performance differences between the generalized and maximum likelihood models in Figure 5 are due to constraining α to $\alpha =1$ (in $P_R\left(y\right|x,\beta )$), or due to inherent difficulty of learning from single counts. The resulting ROC analysis is shown in Figure 6. The data driven setting of hyperparameters $\alpha ,\beta$ leads to slight improvement over $P_{AC}\left(y\right|x)$ and $P_R\left(y\right|x,\beta )$.

Figure 6. ROC curves for test distributions $P_{AC}\left(y\right|x)=P_R\left(y\right|x,\beta \to 0)$ (solid black line) and $P_G\left(y\right|x,{\alpha }_{*},{\beta }_{*})$ (dashed red line). Mean rate of the underlying Poisson source was fixed at $\lambda =5$.

Figure 6. ROC curves for test distributions $P_{AC}\left(y\right|x)=P_R\left(y\right|x,\beta \to 0)$ (solid black line) and $P_G\left(y\right|x,{\alpha }_{*},{\beta }_{*})$ (dashed red line). Mean rate of the underlying Poisson source was fixed at $\lambda =5$.

7. Discussion and Conclusion

Studies of learning algorithms traditionally concentrate on situations where potentially ever increasing number of training examples is available. However, there are situations where only extremely small samples can be used in order to perform an inference. In this contribution we concentrated on extreme case of low count data governed by Poisson distribution, where only a single observation is available. We performed a rigorous theoretical investigation of the appropriateness of various model estimators, based on the single observation. We considered a Bayesian approach along the lines of [2], where the model built on the basis of a single observed count is no longer Poisson, even though we know that the generating source is Poisson (but do not know the mean rate).

We showed that the Bayesian approach is more optimal than the regularized maximum likelihood, in the sense that the expected Kullback–Leibler divergence from the source to the model is smaller for the Bayesian approach. Furthermore, we generalized the original model of [2] to account for possible prior information on expected expression counts. Detailed information theoretic study of learning capabilities of such a generalized model was conducted for the case of low count data. We also quantified the effect of Bayesian averaging on its first two moments.

We demonstrated both theoretically and empirically that the Bayesian model averaging on the generalized model can be potentially beneficial. For large λ, the expected divergence $\mathsf{{\rm Y}}(\lambda ,ϵ)$ of the maximum likelihood estimator from the true Poisson source is dominated by the term

$$\lambda \left(log\lambda -\sum _{x=1}^{\infty }P\left(x\right|\lambda )logx\right)$$

since ${lim}_{\lambda \to \infty }{e}^{-\lambda }(ϵ-\lambda logϵ)=0$. We empirically determined that for $\lambda \ge 10$, $\mathsf{{\rm Y}}(\lambda ,ϵ=1)$ expressed in bits is bounded by $0.7<\mathsf{{\rm Y}}(\lambda ,ϵ=1)<0.8$. Hence, for mean Poisson rates $\lambda \ge 10$, the difference between the expected divergences of the Audic–Claverie and ML estimates from the true source is never less than 0.2 bits and never more than 0.3 bits. In other words,

$$0.2<\mathsf{{\rm Y}}(\lambda ,ϵ=1)-\mathcal{E}\left(\lambda \right)<0.3,\lambda \ge 10$$