# Bayes inference II: random variables

##### Bayes inference II: random variables

E. Canuto, Former faculty, Politecnico di Torino, Torino,Italy

September 7, 2020

Draft

###### Acronyms

RV random variable

PDF   probability density function

PMF probability mass function

MLE  maximum likelihood estimate

MAP maximum a posteriori

LF likelihood function

SD standard deviation

###### Introduction

In part I, the Bayes Theorem was illustrated by several examples under the  assumption of a set $\Omega$ of finite outcomes and their events. We extend the theorem to integer and real random variables (RV) [2]. The integer random variable $K=\left&space;\{&space;0\leq&space;k\leq&space;n&space;\right&space;\}$ was already introduced in Part I to describe the possible outcomes of k heads (H) in n repeated trials of coin tossing. Also the set $\left&space;\{&space;H,T&space;\right&space;\}$ can be associated to the integer RV $X=\left&space;\{&space;0\left&space;(&space;tail&space;\right&space;),1&space;(head)\right&space;\},P\left&space;(&space;X=1&space;\right&space;)=p$.  Consider now a finite sequence $X=\left&space;\{&space;X_{1},...,X_{i},...,X_{n}&space;\right&space;\}$  of independent RV $X_{i}=\left&space;\{&space;0,1&space;\right&space;\}$. The random variable K holds $K=\sum_{i=1}^{n}X_{i}$, and from Part I the likelihood function (LF) of the outcome K=k is the Binomial distribution

$P\left&space;(&space;K=k/p\right&space;)=P\left&space;(&space;K=k;n,p&space;\right&space;)=\begin{pmatrix}n&space;\\&space;k\end{pmatrix}p^{k}\left&space;(&space;1-p&space;\right&space;)^{n-k}\;&space;(1)$

In Part I we assumed that the real number  $0<&space;p<1$ was unknown. Here we assume to know a prior probability density function (PDF) $f\left&space;(&space;p;&space;q_{1},q_{2},...\right&space;)$, which implies that p becomes the outcome of a real RV, denoted with p not to abuse of P; $q_{1},q_{2},...$ denote the hyper-parameters of the PDF. The simplest is the uniform probability density, $f\left&space;(&space;p;&space;q_{1},q_{2},...\right&space;)=1$. Actually, we are looking for a generic PDF with the property that the product of likelihood function and prior PDF (proportional to the posterior PDF, see below) has the same form of the prior to ease computations. The prior PDF is then said to be a conjugate prior of the likelihood function. The conjugate prior of the Binomial distribution is the Beta PDF $f\left&space;(&space;p;&space;\alpha&space;,\beta&space;\right&space;)=\textup{Beta}\left&space;(&space;p;\alpha&space;,\beta&space;\right&space;)$, which is defined for $0\leq&space;p\leq&space;1$ and holds

$\begin{matrix}\textup{Beta}\left&space;(&space;p;\alpha&space;,\beta&space;\right&space;)=\frac{p^{\alpha&space;}\left&space;(&space;1-p&space;\right&space;)^{\beta&space;}}{B\left&space;(&space;\alpha&space;,\beta&space;\right&space;)}&space;\\&space;B\left&space;(&space;\alpha&space;,\beta&space;\right&space;)=\int_{0}^{1}p^{\alpha&space;}\left&space;(&space;1-p&space;\right&space;)^{\beta&space;}dp=\frac{\Gamma&space;\left&space;(&space;\alpha+1&space;\right&space;)\Gamma&space;\left&space;(&space;\beta&space;+1\right&space;)}{\Gamma&space;\left&space;(&space;\alpha+&space;\beta+2&space;\right&space;)}&space;\end{matrix}\:&space;\;&space;(1)$

where $\Gamma&space;\left&space;(&space;k&space;\right&space;)=(k-1)!,&space;\:&space;\;&space;k\:&space;\;&space;\textup{integer}$, -1 has been avoided in the exponents and $B\left&space;(&space;\alpha&space;,\beta&space;\right&space;)$ is the normalizing constant. The uniform density follows from $\alpha&space;=\beta&space;=0$. The PDF in (1) is symmetric around p=0.5 for $\alpha&space;=\beta$. As  the exponents increase, (1) becomes narrower and approaches a Gaussian distribution as Figure 1 shows.


Figure 1 - Symmetric Beta PDF for increasing parameter value.

Remark. Meaning and difference between Binomial and Beta distributions. Comparison of (1) with the Binomial distribution $\begin{matrix}\textup{Bin}\left&space;(&space;k,n,p&space;\right&space;)=P\left&space;(&space;K=k/p,n&space;\right&space;)=\begin{pmatrix}n&space;\\&space;k\end{pmatrix}p^{k}\left&space;(&space;1-p&space;\right&space;)^{n-k}&space;\\&space;\sum_{k=0}^{n}\textup{Bin}\left&space;(&space;k,n,p&space;\right&space;)=1&space;\end{matrix}$

shows their similarity, although they express different models. The binomial RV K is discrete as it accounts for the head counts in n trials, having denoted P(head) (known or unknown) with p. The Beta distribution is the PDF of a bounded real  RV, $0\leq&space;p\leq&space;1$, which we interpret as P(head) .  The Beta exponents may be interpreted as the head and tail counts that has been used to build up the PDF. In fact, as shown in Figure 1, by increasing the exponents, the PDF narrows implying smaller uncertainty.

The goal is to compute the posterior PDF $g\left&space;(&space;p/K=k&space;\right&space;)$ of the parameter p given an outcome k of K obtained from an outcome x of the sequence X. Bayes theorem for random variables requires four probability functions (either discrete or real)

1. the likelihood function of the observed  outcomes, discrete in this case, $P\left&space;(&space;K=k/p\right&space;)$ in (1), conditioned by a generic outcome p of the parameter prior probability,
2. the prior PDF of the real valued parameter p, in (2), depending on known parameters $\left&space;(&space;\alpha&space;,\beta&space;\right&space;)$ called hyper-parameters,
3. the discrete probability $P\left&space;(&space;K=k\right&space;)$ of the observed outcome k whichever be parameter p; it is obtained from the law of total probability (see Part I) extended to discrete and real RV and does not depend on p (it plays the role of a scale factor in Bayes equation)

$P\left&space;(&space;K=k&space;\right&space;)=\int_{0}^{1}P\left&space;(&space;K=k&space;/p\right&space;)f\left&space;(&space;p;\alpha&space;,\beta&space;\right&space;)dp\:&space;\;&space;(3)$The posterior PDF holds

$\begin{matrix}g\left&space;(&space;p/K=k&space;\right&space;)=f\left&space;(&space;p;\alpha&space;,\beta&space;\right&space;)\frac{P\left&space;(&space;K=k/p\right&space;)}{\int_{0}^{1}P\left&space;(&space;K=k&space;/p\right&space;)f\left&space;(&space;p;\alpha&space;,\beta&space;\right&space;)dp}&space;\\&space;=\frac{p^{\alpha+k&space;}\left&space;(&space;1-p&space;\right&space;)^{\beta+n-k&space;}}{B\left&space;(&space;\alpha&space;+k,\beta&space;+n-k\right&space;)}\end{matrix}\;&space;(4)$

In Part I , the likelihood function was maximized to provide the MLE of p equal to $\hat{p}=k/n$. The same can be done here, by maximizing the posterior probability versus the unknown p. The new estimate, known as the Maximum a Posteriori (MAP) estimate holds

$\begin{matrix}\frac{d}{dp}g\left&space;(&space;p/K=k&space;\right&space;)=p^{\alpha&space;+k-1}\left&space;(&space;1-p&space;\right&space;)^{\beta&space;+n-k-1}\left&space;(\left&space;(\alpha&space;+k&space;\right&space;)&space;\left&space;(&space;1-p&space;\right&space;)-\left&space;(&space;\beta&space;+n-k&space;\right&space;)p&space;\right&space;)\\&space;\hat{p}_{MAP}=\frac{\alpha&space;+k}{\alpha&space;+\beta&space;+n}\end{matrix}\;&space;(5)$

The MAP estimate becomes the MLE for $\alpha&space;=\beta&space;=0$, namely under uniform prior density, which is intuitive as the uniform density does not depend on p. We can also say that (5) is the mode of the posterior PDF.  The posterior Beta distribution in (4) for a prior with $\alpha&space;=\beta&space;=1$ and data n=20 and k={5,10} is shown in Figure 2. The posterior PDF becomes narrower than prior, meaning that uncertainty is reduced. Although the prior is centered in p=0.5, the posterior becomes centered below p=0.5 when k=5 (biased coin), a reasonable result.

Figure2 - Prior and two posterior Beta densities.



Remark 1. Posterior mean and variance.  From (4), the posterior PDF is the function $\textup{Beta}\left&space;(&space;\alpha&space;+k,\beta&space;+n-k&space;\right&space;)$, whose mean and variance hold

$\begin{matrix}E\left&space;\{&space;p/K=k&space;\right&space;\}=\frac{\alpha&space;+k+1}{\alpha&space;+\beta&space;+n+2},&space;\:&space;\;&space;\lim_{n,k\rightarrow&space;\infty&space;}E\left&space;\{&space;p/K=k&space;\right&space;\}&space;=k/n\\&space;\lim_{n\rightarrow&space;\infty&space;}&space;var\left&space;\{&space;p/K=k&space;\right&space;\}&space;=\frac{k\left&space;(&space;n-k&space;\right&space;)}{n^{3}}=\frac{\hat{p}_{MLE}\left&space;(&space;1-\hat{p}_{MLE}&space;\right&space;)}{n}&space;\end{matrix}$

The limiting posterior variance depends on the MLE estimation and the size n of the repetitions (size of the measurements), implying that the estimate uncertainty goes down with $1/\sqrt{n}$.

Remark 2. Iterative posterior. Posterior PDF and MAP estimate can be easily improved when new data (new coin tossing outcomes, $(\Delta&space;k,\Delta&space;n)$) become available. It is enough to replace $(k,n)$ with $(k+\Delta&space;k,n+\Delta&space;n)$.

Remark 3. Prediction of a future outcome and zero count paradox. Given a posterior PDF of $p=P\left&space;(&space;\textup{head}\right&space;)=P\left&space;(&space;X_{i}=1&space;\right&space;)$,

$g\left&space;(&space;p/K=k&space;\right&space;)=\textup{Beta}\left&space;(&space;a,b&space;\right&space;);&space;a=\alpha&space;+k,&space;b=\beta&space;+n-k$we want the probability $P\left&space;(&space;X_{n+1}=1/K=k&space;\right&space;)$ that the next trial n+1 provides head as the outcome. To this end, we employ the law of total probability

$\begin{matrix}P\left&space;(&space;X_{n+1}=1/K=k&space;\right&space;)=\int_{0}^{1}P(X_{n+1}=1/p)g\left&space;(&space;p/K=k&space;\right&space;)dp=&space;\\&space;=\int_{0}^{1}p\textup{Beta}\left&space;(&space;\alpha&space;+k,\beta&space;+n-k&space;\right&space;)dp=E\left&space;\{&space;p/K=k&space;\right&space;\}=\frac{\alpha&space;+k+1}{\alpha&space;+\beta&space;+n+2}&space;\end{matrix}\:&space;(6)$

The result in (6) is equivalent to the identity $P\left&space;(&space;X_{n+1}=1&space;\right)=E\left&space;\{&space;p/K=k&space;\right&space;\}$ in (6).   Alternatively, one may adopt the MLE instead of the posterior mean, that is: $P\left&space;(&space;X_{n+1}=1&space;\right)=\hat{p}_{MLE}\left&space;(&space;k,n&space;\right&space;)=k/n$. Which is the difference?  Consider the case when k=0, which may occur under small n. The MLE fails to provide a non zero estimate of p, a paradox known as the zero count problem, where inference (induction) from observations fails. On the contrary, the posterior mean provides a non zero estimate of p, also when the prior probability of p is uniform, namely $\alpha&space;=\beta&space;=0$, in which case from (6) we have

$\begin{matrix}\hat{p}_{MAP}\left&space;(&space;\alpha&space;=\beta&space;=k=0&space;\right&space;)=\frac{1}{n+2}&space;\\&space;\hat{p}_{MAP}\left&space;(&space;\alpha&space;=\beta&space;=0&space;\right&space;)=\frac{k+1}{n+2}&space;\;&space;(\textup{Laplace's&space;rule&space;of&space;succession})&space;\end{matrix}$

Let us remark that under uniform prior and k=0, also the MAP estimate in (5) would fail to provide a non zero estimate, which fact suggests the practice, known as add-one smoothing, of adding 1 to empirical counts. Addition of 2 in the denominator of (6) is justified by adding 1 to both head count k and tail count n-k.

Remark 4. Prediction of future repeated trials

We aim to predict the number h of heads in m future repeated trials, after having observed k heads in n trials and having assumed a prior PDF of p. We define the random variable $H=\sum_{i=n+1}^{n+m}X_{i}$. We repeat (6) by replacing   $P(X_{n+1}=1/p)$ with the binomial distribution with parameters $(m,h)$

$\begin{matrix}P\left&space;(&space;H=h/K=k&space;\right&space;)=\int_{0}^{1}P(H=h/p)g\left&space;(&space;p/K=k&space;\right&space;)dp=&space;\\&space;=\begin{pmatrix}&space;m&space;\\&space;h\end{pmatrix}&space;\int_{0}^{1}p^{h}\left&space;(&space;1-p&space;\right&space;)^{m-h}\textup{Beta}\left&space;(&space;\alpha&space;+k,\beta&space;+n-k&space;\right&space;)dp=\\=&space;\begin{pmatrix}&space;m&space;\\&space;h\end{pmatrix}\frac{1}{B(\alpha&space;+k,\beta&space;+n-k&space;)}&space;\int_{0}^{1}p^{h+\alpha&space;+k}\left&space;(&space;1-p&space;\right&space;)^{m-h+\beta&space;+n-k}dp&space;\\&space;=\begin{pmatrix}&space;m&space;\\&space;h\end{pmatrix}\frac{B(h+\alpha&space;+k,m-h+\beta&space;+n-k&space;)}{B(\alpha&space;+k,\beta&space;+n-k&space;)}&space;\end{matrix}\:&space;(7)$The last expression as a function of h is an integer probability function known as beta-binomial, $\textup{BetaBin}(h/m,n,k,\alpha&space;,\beta&space;))$. Of course  $\sum_{h=0}^{m}P\left&space;(&space;H=h/K=k&space;\right&space;)=1$. The mean value holds

$E\left&space;\{&space;H=h&space;\right&space;\}=m\frac{\alpha&space;+k+1}{\alpha&space;+\beta&space;+n+2}$

which for m=1 provides the same expected value as in (6). Figure 3 shows the prediction probability of head count during future m=10 trials under two batches of observed data: (k=9,n=20) and (k=45,n=100). The prior parameters are $\alpha&space;=\beta&space;=1$.  In the second case, the prediction probability becomes narrower.

Figure 3 - Prediction probability of head counts out of m=10 trials, after n=(20,100) observed trials with head count k=(8,45).

Application to ballot

Head and tail outcomes can be replaced by ballot candidates (C1,C2). Let us keep the RV notation $X_{i}=\left&space;\{&space;T=0&space;(\textup{C1}),&space;H=1&space;(\textup{C2}))&space;\right&space;\}$.  A repeated trial of size n can be either an election poll or the same election. Election polls are the available observations, which provide candidate preference counts (k,n-k). We aim to compute the posterior probability that a candidate will be elected, given the poll data.  Let us denote with p the probability that a vote is for candidate 2, $P\left&space;(&space;X_{i}&space;=1\right&space;)=p$. The probability that candidate 2 will win is $P\left&space;(&space;p>1/2&space;\right&space;)$, which implies that p is a RV, as in the Bayes inference. As an alternative, we can get the MLE estimate $p\hat{}_{MLE}=k/n$, which would imply the election of C2 under $p\hat{}_{MLE}>1/2$. Unfortunately the estimate is uncertain (it is a RV) and it must be completed  by its variance. Based on data in [1] we will employ both methods.  The difficulty of Bayesian inference  is the prior PDF of p.  One way  is to subdivide the available poll data in two parts: the oldest will estimate the prior PDF, the newest will estimate the likelihood function.  Poll data form the same company are shown in Table 1.

 Table I – Election polls, $i=1,...,4$ No Date Candidate 1, $n_{i}-k_{i}$,$X_{i}=0$ Candidate 2, $k_{i}$, $X_{i}=1$ Total, $n_{i}$ Used for 1 September 4-7 344 284 628 Prior PDF 2 September 25-28 325 312 637 Prior PDF 3 October 17-20 339 346 685 Prior PDF 4 Partial sum $\beta+1=1008$ $\alpha+1=942$ 1950 Prior PDF 5 Last poll $n_{4}-k_{4}=511$ $k_{4}=556$ 1067 Likelihood 6 November 2 election (fraction) 0.511 0.489 1

The exponents $\left&space;(\alpha+1=\sum_{i=1}^{3}k_{i}&space;,\beta+1=\sum_{i=1}^{3}\left&space;(n_{i}-k_{i}&space;\right&space;)&space;\right&space;)$ of the prior  PDF are estimated from the sum of the preferences in the oldest three polls (row 4 in table 1). The exponents of the likelihood function correspond to the last poll  preferences. The posterior PDF in (4) is a Beta PDF with parameters  $\left&space;(\alpha+k_{4}&space;,\beta+n_{4}-k_{4}&space;\right&space;)$. The prior PDF (blue line), the normalized LF (normalized to become a Beta PDF, red line) and the posterior PDF (green line) are in Figure 4. Due to large batch of data (about 1000 people), the three functions being very narrow, have been plotted for $0.4\leq&space;p\leq&space;0.6$. The prior PDF is in favor of C1, the LF from the last poll is in favor of C2, but the posterior PDF is in favor of C1, but to a less extent than the prior (see Table 2, below). The posterior SD is smaller than the prior. The posterior mean value in Remark 1 and the standard deviation (SD) hold

$E\left&space;\{&space;g\left&space;(&space;p/K=k&space;\right&space;)&space;\right&space;\}\cong&space;0.4965,&space;SD\left&space;\{&space;g&space;\left&space;(&space;p/K=k&space;\right&space;)\right&space;\}\cong&space;0.0091$

The posterior probability that  $p>1/2$ (C2 probability of being voted) holds $P\left&space;(&space;p>1/2&space;\right&space;)\cong&space;0.351$, a rather small value due to mass concentration close to mean value, which is coherent  with the SD.  This implies that, though LF looks favorable to C2, the posterior PDF turns in favor of C1, because of the unfavorable prior.  However, as proved below, the MLE computed from four poll data perfectly agrees with Bayesian inference results, as expected, due to large data batches, and the fact that the posterior PDF is proportional to the four poll MLE and their mean value and mode coincide.


Figure 4 - The components of the Bayesian inference: the prior Beta PDF (in blue, from the oldest polls), the normalized likelihood function  (in red, from the last poll) and the posterior Beta whose mode is less than 1/2.                        

A similar result can be found from the cumulative MLE of the four poll data in Table 1:

$\hat{p}_{\textup{MLE,&space;cum}}=\frac{\sum_{i=1}^{4}k_{i}}{\sum_{i=1}^{4}n_{i}}\cong&space;0.4965<0.5$

which is the same value of the posterior mean. Figure 5 shows the cumulative MLE (horizontal green line) together with the sequence of the single poll MLEs. The winning probability of the candidate 2 (C2) is increasing with the ballot campaign, but the cumulative MLE shows that they look insufficient to reach a winning probability larger than 1/2, a result in agreement with Bayesian inference. Data from [1] refer to 2004 US presidential elections, but restricted to Ohio state (November 2004). Actually they were not a ballot, but the score of other candidates was negligible. C2 did not win the ballot (see row 6 of Table 1).

Figure 5 - History of the winning probability of the candidate 2 (C2) from single poll data (pointwise blue) and cumulative data (pointwise green). The 1/2 threshold is in red.

Bayesian inference and MLE results are summarized in Table 2.

 Table 2 – Prior, posterior poll estimates of C2 preference probability No Parameter Symbol Unit Value Remark 1 Prior mean $E\left&space;\{&space;f\left&space;(&space;\right&space;)&space;\right&space;\}$ fraction 0.4830 From i=1,2,3 2 Prior SD $SD\left&space;\{&space;f\left&space;(&space;\right&space;)&space;\right&space;\}$ fraction 0.0113 same 3 Posterior mean $E\left&space;\{&space;g\left&space;(&space;\right&space;)&space;\right&space;\}$ fraction 0.4965 From i=1,2,3,4 4 Posterior SD $SD\left&space;\{&space;g\left&space;(&space;\right&space;)&space;\right&space;\}$ fraction 0.0091 same 5 Last poll MLE $\hat{p}_{MLE,4}$ fraction 0.521 From i=4 6 Cumulative poll MLE $\hat{p}_{MLE,cum}$ fraction 0.4965 From i=1,2,3,4 7 Voting result $p_{true}$ fraction 0.489

Remark. The progression of the MLE estimates in Figure 5 in contrast with the voting result (row 7 of Table 2) looks a bit odd, since the latter appears more coherent with the oldest polls than those close to election day. The confirmation comes from the cumulative poll MLE (row 6, Table 2) which shows coherency with the voting result, but at the price of a uniform weight of past and recent data. We cannot investigate further this oddity, but clearly poll data hide other assumptions than those employed here.

###### Multinomial inference

In the previous section Bayesian inference aimed to infer the occurrence probability of one outcome between two alternatives (binomial inference) like {head, tail}, {candidate 1, candidate 2}. Here we extend the inference to an arbitrary finite set $\left&space;\{&space;\omega&space;_{1},...,\omega&space;_{k},&space;...,\omega&space;_{m}&space;\right&space;\}$ of outcomes, like for instance the m=6 faces of a die. Let us associate the integer RV $X_{i}=\left&space;\{&space;1,&space;..,k,...,m&space;\right&space;\}$ with generic outcome $x_{i}$. As in the binomial case, we assume to observe n independent repetitions of the same experiment, or which is the equivalent, n extractions $X=&space;\left&space;\{&space;X_{i},i=1,...,n&space;\right&space;\}$ of $X_{i}$, whose generic outcome is denoted by $x=\left&space;\{&space;x_{1},...,x_{i},...,x_{n}&space;\right&space;\}$. Let us denote with $p_{k}$ the occurrence probability of the outcome k, such that

$\sum_{k=1}^{m}p_{k}=1;&space;0\leq&space;p_{k}\leq&space;1\;&space;(8)$

and with $n_{k}$  the number of  occurrences in the observed sequence x.  The likelihood function of the repeated trials is the probability of observing $n_{k},&space;k=1,...,m$ repetitions of the outcome k in a specific sequence x. Let us introduce the inline vector $\mathbf{p}=\left&space;[&space;p_{1},...,p_{k},...,p_{m}&space;\right&space;]$. It holds

$P\left&space;(&space;X=x/\mathbf{p}&space;\right&space;)=\prod_{k=1}^{m}p_{k}^{n_{k}}\:&space;\;&space;(9)$

Remark 1. MLE. From (9) we can compute the MLE, by taking the logarithm of the LF and adding the constraint in (8) multiplied by a Lagrange multiplier $\lambda$, as follows

$\begin{matrix}L(&space;\mathbf{p},\lambda&space;)=\textup{log}P\left&space;(&space;X=x/\mathbf{p}&space;\right)+\lambda&space;\left&space;(&space;1-\sum_{k=1}^{m}&space;p_{k}\right&space;)=&space;\\&space;=\sum_{k=1}^{m}n_{k}\textup{log}p_{k}+\lambda&space;\left&space;(&space;1-\sum_{k=1}^{m}&space;p_{k}\right&space;)&space;\end{matrix}\:&space;\;&space;(10)$

Setting the derivatives of the Lagrangian in (10) to zero we obtain

$\begin{matrix}\frac{\partial&space;L}{\partial&space;p_{k}}=\frac{n_{k}}{p_{k}}-\lambda=0;&space;\;&space;\frac{\partial&space;L}{\partial&space;\lambda&space;}=1-\sum_{k=1}^{k}p_{k}=0&space;\\&space;\frac{n_{k}}{p_{k}}=\lambda&space;\Rightarrow&space;n=\sum_{k=1}^{k}n_{k}=\lambda\sum_{k=1}^{k}p_{k}&space;=\lambda&space;\\&space;\hat{p}_{k,MLE}=\frac{n_{k}}{\lambda&space;}=\frac{n_{k}}{n&space;}\geqslant&space;0&space;\end{matrix}\:&space;\;&space;(10)$

namely the MLE which satisfies the constraints in (8).

The conjugate prior of the LF is the multinomial Beta PDF, known as Dirichlet PDF (see the Appendix). By avoiding -1 in the exponent, by defining the inline vector $\mathbf{p}=\left&space;[&space;p_{1},...,p_{k},...,p_{m}\right&space;]$  and the exponent vector $\mathbf{a&space;}=\left&space;[&space;\alpha&space;_{1},..,\alpha&space;_{k},&space;...,\alpha&space;_{m}\right&space;]$, we write

$\begin{matrix}\textup{Dir}(\mathbf{p};\mathbf{a})=f\left&space;(\mathbf{&space;p};\mathbf{a}&space;\right&space;)=\frac{1}{B\left&space;(\mathbf{a}&space;\right&space;)}\prod_{k=1}^{m}p_{k}^{\alpha&space;_{k}}&space;\\B\left&space;(\mathbf{a}&space;\right&space;)=\frac{\prod_{k=1}^{m}\Gamma&space;\left&space;(\alpha&space;_{k}&space;+1\right&space;)}{\Gamma&space;\left&space;(&space;\sum_{k=1}^{m}&space;\left&space;(\alpha_{k}+1&space;\right&space;)\right&space;)}&space;\end{matrix}\;&space;\left&space;(11&space;\right&space;)$

Composition of prior PDF and likelihood function provides, as in the Introduction, the posterior PDF in the Dirichlet form as follows

$\begin{matrix}\textup{Dir}(\mathbf{p};\mathbf{a+n})=g\left&space;(\mathbf{&space;p}/X=x;\mathbf{a+n}&space;\right&space;)=\frac{1}{B\left&space;(\mathbf{a+n}&space;\right&space;)}\prod_{k=1}^{m}p_{k}^{\alpha&space;_{k}+n_{k}}&space;\\B\left&space;(\mathbf{a+n}&space;\right&space;)=\frac{\prod_{k=1}^{m}\Gamma&space;\left&space;(\alpha&space;_{k}+n_{k}&space;+1\right&space;)}{\Gamma&space;\left&space;(&space;\sum_{k=1}^{m}&space;\left&space;(\alpha&space;_{k}+n_{k}+1&space;\right&space;)\right&space;)}&space;\end{matrix}\;&space;\left&space;(12&space;\right&space;)$

Remark 2. MAP estimate. It is obtained in the same way as the MLE by maximizing the logarithm of (11) added with the constraint in (8) multiplied by $\lambda$. The estimate follows from (10) by replacing $n_{k}$ with $n_{k}+\alpha&space;_{k}$, that is

$\hat{p}_{k,MAP}=\frac{n_{k}+\alpha&space;_{k}}{n+\sum_{k=1}^{m}\alpha&space;_{k}&space;}&space;\:&space;\;&space;(13)$

Remark. Zero count problem. The MLE  estimate provides zero estimate when $n_{k}=0$. The MAP estimate provides zero estimate when in the prior PDF is uniform and  $\alpha&space;_{k}=0$.

###### Naive Bayesian classifiers

Consider n discrete independent RV $\left&space;\{&space;X_{k},k=1,...,m&space;\right&space;\}$, called features which are employed to classify into n classes/categories $C=\left&space;\{&space;c_{i}&space;,&space;i=1,.,n\right&space;\}$ elements of a population, like documents, people, generic objects, … Let us denote an outcome of the feature RV with $\mathbf{x}=[&space;x_{1},...,x_{k},...,x_{m}&space;]$. Each feature may assume a finite, discrete but arbitrary set of values, $X_{k}=\left&space;\{&space;x_{k1},&space;...,x_{kj},...,x_{kn(k)}&space;\right&space;\}$.  Given the Cartesian product $X=X_{1}\times&space;...\times&space;X_{k}\times&space;...\times&space;X_{m}$, a class $c_{i}\subset&space;X,\;&space;c_{i}\cap&space;c_{j}=0,i\neq&space;j$ is a known subset of  $X$, such that $\bigcup_{i=1}^{n}c_{i}=X$. More generically features may also assume real values. We aim to find the posterior probability $P\left&space;(&space;C=c_{i}&space;/X=\mathbf{x}\right&space;)$ of classifying the feature outcome $\mathbf{x}$ into $c_{i}$. The posterior can be rewritten by the Bayes theorem as

$P\left&space;(&space;C=c_{i}&space;/X=\mathbf{x}\right&space;)=P\left&space;(&space;C=c_{i}&space;\right&space;)\frac{P\left&space;(&space;X=\mathbf{x}/C=c_{i}&space;\right&space;)}{P\left&space;(&space;X=\mathbf{x}&space;\right&space;)}=\frac{P\left&space;(&space;C=c_{i},X=\mathbf{x}&space;\right&space;)}{P\left&space;(&space;X=\mathbf{x}&space;\right&space;)}$

where $P\left&space;(&space;C=c_{i},X=\mathbf{x}&space;\right&space;)$ is the joint probability that the feature outcome $\mathbf{x}$ and the class $c_{i}$ occurs together. The naive assumption amounts to how compute this joint probability. Let us factor the joint probability by repeated used of conditional probabilities

$\begin{matrix}P\left&space;(&space;C=c_{i},X=\mathbf{x}&space;\right&space;)=P\left&space;(&space;x_{1}/C=c_{i},x_{2},...,x_{m}\right&space;)P\left&space;(&space;C=c_{i},x_{2},...,x_{m}\right&space;)=&space;\\&space;=P\left&space;(&space;x_{1}/C=c_{i},x_{2},...,x_{m}\right&space;)P\left&space;(&space;x_{2}/C=c_{i},x_{3},...,x_{m}\right&space;)P\left&space;(&space;C=c_{i},x_{3},...,x_{m}\right&space;)=&space;\\&space;...&space;\\&space;=P\left&space;(&space;x_{1}/C=c_{i},x_{2},...,x_{m}\right&space;)...P\left&space;(&space;x_{m-1}/C=c_{i},x_{m}\right&space;)P\left&space;(&space;x_{m}/C=c_{i}\right&space;)P\left&space;(&space;C=c_{i}\right&space;)&space;\end{matrix}$

The naive assumption simplifies the above expression, by assuming that given a category $c_{i}$, the m different features occur independently one from another (conditional independence). In other words, when we observe features related to different categories they may be correlated, but when we observe features within a well defined category, they are independent. In general the assumption is far form reality, except when the values of class features are variations around a mean value, typical of the class, and variations are feature uncorrelated.

For instance, you want to classify a day from temperature, wind, humidity, outlook (sunny, cloudy, overcast, rainy, …) in view of open air activities. Classes may be: good, uncertain, bad. Of course the above features are in general correlated, but we assume them independent within a class.

###### Reference

[1] Basics of Bayesian Statistics, from http://www.stat.cmu.edu/~brian/463-663/week09/Chapter%2003.pdf.

[2] K.P. Murphy, Machine learning: a probabilistic perspective, The MIT Press, Cambridge, MA, 2012.

###### Appendix: discrete PMF and real PDF

We will use the Iverson bracket, which generalizes the Kronecker delta. Given a statement S which may be either true or false, we write

$\left&space;[&space;S&space;\right&space;]=\left\{\begin{matrix}1,&space;\textup{true}&space;\\&space;0,&space;\textup{false}&space;\end{matrix}\right.$

Discrete  probability mass functions(or distributions)

Bernoulli distribution $\textup{Ber}(x;p)$. Given a binary RV $X=\left&space;\{&space;0,1&space;\right&space;\}$ , the generic outcome x and the parameter $0\leq&space;p\leq&space;1$, the PMF is defined as

$\begin{matrix}\textup{Ber}(x;p)=P\left&space;(&space;X=x&space;\right&space;)=p^{\left&space;[&space;x=1&space;\right&space;]}\left&space;(&space;1-p&space;\right&space;)^{\left&space;[&space;x=0&space;\right&space;]}&space;\\&space;E\left&space;\{&space;X&space;\right&space;\}=p,\;&space;\textup{var}\left&space;\{&space;X&space;\right&space;\}=p\left&space;(&space;1-p&space;\right&space;)&space;\end{matrix}&space;\:&space;\;&space;(A1)$

Binomial distribution $\textup{Bin}(x;p,k,n)$. Given a sequence $X=\left&space;\{&space;X_{i},i=1...,n&space;\right&space;\}$ of independent binary RV $X_{i}=\left&space;\{&space;0,1&space;\right&space;\}$  we construct the RV $K=\sum_{i=1}^{n}X_{i}$ with generic outcome $0\leq&space;k\leq&space;n$. The PMF is defined as

$\begin{matrix}\textup{Bin}(k;n,p)=P\left&space;(&space;K=k&space;\right&space;)=\begin{pmatrix}n\\&space;k\end{pmatrix}p^{k}\left&space;(&space;1-p&space;\right&space;)^{n-k}&space;\\&space;E\left&space;\{&space;K&space;\right&space;\}=np,\;&space;\textup{var}\left&space;\{&space;K&space;\right&space;\}=np\left&space;(&space;1-p&space;\right&space;)&space;\end{matrix}\;&space;(A2)$

Categorical distribution $\textup{Cat}(x;\mathbf{p},m),&space;\mathbf{p}=[p_{1},...,p_{k},...,p_{m}]$Given a discrete RV $X=\left&space;\{&space;x_{1},...,x_{k},...,x_{m}&space;\right&space;\}$ where the integer $x_{k}$ denotes a class or category of different objects or features (e.g., die faces) with probability $0\leq&space;p_{k}\leq&space;1$ such that $\sum_{k=1}^{m}p_{k}=1$, the PMF is defined by

$\begin{matrix}\textup{Cat}(x;\mathbf{p})=P\left&space;(&space;X=x_{k}&space;\right&space;)=\prod_{k=1}^{m}p_{k}^{\left&space;[&space;x=x_{k}&space;\right&space;]}=\sum_{k=1}^{m}&space;\left&space;[&space;x=x_{k}&space;\right&space;]p_{k}\\&space;\end{matrix}\;&space;\left&space;(A3&space;\right&space;)$

Multinomial distribution $\textup{Mu}(x;n,m,&space;\mathbf{p}),&space;\mathbf{p}=[p_{1},...,p_{k},...,p_{m}]$. Given a sequence of  $X=\left&space;\{&space;X_{i},i=1...,n&space;\right&space;\}$ of n identical and independent categorical RV  $X_{i}=\left&space;\{&space;x_{1},...,x_{k},...,x_{m}&space;\right&space;\}$, the PMF of the random vector $\mathbf{\mathrm{}N}=\left&space;[&space;N_{1},...,N_{k},...,N_{m}\right&space;]$ ,where the RV $N_{k}$ accounts for the occurrence number of the feature $x_{k}$,

$\begin{matrix}\textup{Mu}(\mathbf{n};n,m,\mathbf{p})=P\left&space;(&space;\mathbf{N}=\mathbf{n}&space;\right&space;)=\begin{pmatrix}n\\&space;n_{1}...n_{k}...n_{m}\end{pmatrix}\prod_{k=1}^{m}p_{k}^{n_{k}}&space;\\&space;E\left&space;\{&space;N_{}&space;\right&space;\}=np_{k},\;&space;\textup{var}\left&space;\{&space;N_{k}&space;\right&space;\}=np_{k}\left&space;(&space;1-p_{k}&space;\right&space;)&space;\end{matrix}\:&space;\;&space;(A4)$

The categorical PMF is equivalent to the multinomial distribution with n=1, that is

###### $\textup{Cat}(x;\mathbf{p})=\textup{Mu}(\mathbf{n};1,m,\mathbf{p});&space;\;&space;\mathbf{n}=[0&space;...1...0]\:&space;\;&space;(A5)$

Probability density functions

Gaussian or normal distribution. The univariate normal PDF holds

$\begin{matrix}N\left&space;(&space;x;\mu&space;,\sigma&space;\right&space;)=\frac{1}{\sigma&space;\sqrt{2\pi&space;}}\textup{exp}\left&space;(&space;-\frac{\left&space;(&space;x-\mu&space;\right&space;)^{2}}{2\sigma&space;^{2}}&space;\right&space;)&space;\\&space;E\left&space;\{&space;x&space;\right&space;\}=\mu&space;,\textup{var}\left&space;\{&space;x&space;\right&space;\}=\sigma&space;^{2}&space;\end{matrix}\:&space;\;&space;(A6)$

The standard normal PDF assumes $\mu&space;=0,\sigma&space;=1$. The multivariate Gaussian PDF of a RV $\mathbf{x}=\left&space;[&space;x_{1},...,x_{i},...,x_{n}&space;\right&space;]$ holds

$\begin{matrix}N\left&space;(&space;\mathbf{x};\mathbf{\mu}&space;,S^{2}&space;\right&space;)=\frac{1}{\left&space;(2\pi&space;\right&space;)^{n/2}\textup{det}S}\textup{exp}\left&space;(&space;-\frac{1}{2}\left&space;(&space;\mathbf{x}-\mathbf{\mu}&space;\right&space;)^{T}S&space;^{-2}&space;\left&space;(&space;\mathbf{x}-\mathbf{\mu}&space;\right&space;)\right)&space;\\&space;E\left&space;\{&space;\mathbf{x}&space;\right&space;\}=\mathbf{\mu}&space;,\textup{covar}\left&space;\{&space;\mathbf{x}&space;\right&space;\}=S^{2}&space;\end{matrix}\;&space;(A7)$

where the covariance matrix has been raised to the exponent 2 being a nonnegative definite matrix, with nonnegative eigenvalues and (not unique) square root such that $S^{2}=S^{T}S$. $S&space;^{2}\left&space;(&space;i,i&space;\right&space;)$ is the variance of $x_{i}$, and $\textup{tr}S^{2}=\sum_{i=1}^{n}S^{2}(i,i)$ equals the sum of the eigenvalues.

Beta distribution. The univariate Beta PDF refers to the bounded RV $0\leq&space;X\leq&space;1$. By dropping -1 from the exponents, we write

$\begin{matrix}\textup{Beta}\left&space;(&space;x;\alpha,&space;\beta&space;\right&space;)=\frac{x^{\alpha&space;}\left&space;(&space;1-x&space;\right&space;)^{\beta&space;}}{B\left&space;(&space;\alpha&space;,\beta&space;\right&space;)},\;&space;B\left&space;(&space;\alpha&space;,\beta&space;\right&space;)=\frac{\Gamma&space;\left&space;(&space;\alpha&space;+1&space;\right&space;)\Gamma&space;\left&space;(&space;\beta&space;+1&space;\right&space;)}{\Gamma&space;\left&space;(&space;\alpha&space;+\beta+2&space;\right&space;)}&space;\\&space;E\left&space;\{&space;x&space;\right&space;\}=&space;\frac{\alpha&space;+1}{\alpha+\beta&space;+2};&space;\;\textup{&space;var}\left&space;(&space;x&space;\right&space;)=\frac{\left&space;(&space;\alpha&space;+1&space;\right&space;)\left&space;(&space;\beta&space;+1&space;\right&space;)}{\left&space;(&space;\alpha&space;+\beta+2&space;\right&space;)^{2}\left&space;(&space;\alpha&space;+\beta+3&space;\right&space;)}&space;\end{matrix}\:&space;\;&space;(A8)$

The multivariate Beta PDF, known as Dirichlet PDF, is defined over the simplex $S_{m}=\left&space;\{&space;\mathbf{x};0&space;\leq&space;x_{k}\leq&space;1,\sum_{k=1}^{m}&space;x_{k}=1\right&space;\}$. The PDF expression is

$\begin{matrix}\textup{Dir}\left&space;(&space;\mathbf{x};&space;\mathbf{a&space;}=\left&space;[&space;\alpha&space;_{1},...,\alpha&space;_{m}&space;\right&space;]\right&space;)=\frac{\prod_{k=1}^{m}x_{k}^{\alpha&space;_{k}}\left&space;[&space;\mathbf{x}\in&space;S_{m}&space;\right&space;]}{B\left&space;(&space;\mathbf{a}&space;\right&space;)},\;&space;B\left&space;(&space;\mathbf{a}&space;\right&space;)=\frac{\prod_{i=1}^{m}\Gamma&space;\left&space;(&space;\alpha&space;_{i}+1\right&space;)}{\Gamma&space;\left&space;(&space;\sum_{i=1}^{m}\left&space;(&space;\alpha&space;_{k}+1\right&space;)&space;\right&space;)}&space;\\&space;E\left&space;\{&space;x_{k}&space;\right&space;\}=\frac{\alpha&space;_{k}+1}{\alpha_{0}&space;+m&space;},\;&space;\textup{var}\left&space;\{&space;x_{k}&space;\right&space;\}=\frac{\left&space;(\alpha&space;_{k}+1&space;\right&space;)\left&space;(&space;\alpha_{0}&space;+m+&space;1-\alpha&space;_{k}&space;\right&space;)}{\left&space;(&space;\alpha_{0}&space;+m&space;\right&space;)^{2}\left&space;(&space;\alpha_{0}&space;+m+1&space;\right&space;)&space;},\:&space;\;&space;\alpha_{0}&space;=\sum_{i=1}^{m}\alpha&space;_{k}&space;\end{matrix}\;&space;(A9)$

Under uniform exponent $\alpha&space;_{k}=\alpha&space;/m$, mean and variance in (A9) become

$E\left&space;\{&space;x_{k}&space;\right&space;\}=\frac{1}{m&space;},\;&space;\textup{var}\left&space;\{&space;x_{k}&space;\right&space;\}=\frac{1}{m^{2}}\frac{1-1/m}{&space;\alpha&space;+1+1/m&space;}&space;\;&space;(A10)$

Thus increasing $\alpha$ and m increases precision.

Remark A1. Let us denote the random vector of the Dirichlet PDF with $X=\left&space;[&space;X_{1},...,X_{k},...,X_{m}&space;\right&space;]$$\mathbf{x}$ denotes the vector of outcomes.   The vector is constrained by $\sum_{k=1}^{m}X_{k}=1$.  Under m=2, $\textup{Dir}\left&space;(&space;\mathbf{x=}\left&space;[&space;x_{1},x_{2}=1-x_{1}&space;\right&space;];&space;\mathbf{a&space;}=\left&space;[&space;\alpha&space;_{1},\alpha&space;_{2}&space;\right&space;]\right)$ is a Beta distribution $\textup{Beta}\left&space;(&space;x_{1}=1-x_{2};\alpha_{1},&space;\alpha_{2}&space;\right&space;)$.

Remark 2. Stick-breaking approach to Dirichlet RV generation. Consider a segment of length y=1 and randomly divided into two intervals of length $x_{1}=y_{1},&space;1-y_{1}$. Then further subdivide the segment of length $1-y_{1}$  into two segments of fractional length $y_{2},1-y_{2}$, which implies that segment k=2 has length $x_{2}=y_{2}\left&space;(&space;1-y_{1}&space;\right&space;)$ and the remaining segment has length $r_{2}=\left&space;(&space;1-y_{2}&space;\right&space;)\left&space;(&space;1-y_{1}&space;\right&space;)$. In other terms, $y_{2}$ is the length of the partition of  a unit segment. By proceeding in the same way, the last segment k=m-1 to be further subdivided, has length $r_{m-2}=\prod_{k=1}^{m-2}\left&space;(&space;1-y_{k}&space;\right&space;)$ . The length of the last two segments holds

$x_{m-1}=y_{m-1}\prod_{k=1}^{m-2}&space;\left&space;(&space;1-y_{k}&space;\right&space;),\;&space;r_{m-1}=\prod_{k=1}^{m-1}&space;\left&space;(&space;1-y_{k}&space;\right&space;)$

The length sum of the m segments holds one, namely $\sum_{k=1}^{m-1}x_{k}+r_{m-1}=1$. Since the sequence $\left&space;\{&space;y_{k}&space;,k=1,...,m-1\right&space;\}$ has been obtained by partitioning segments of unit length, it can be obtained as the outcome of a sequence of Beta  distributions

$\textup{Beta}\left&space;(&space;y_{1};\alpha_{1},\sum_{&space;k=2}^{m}\alpha&space;_{k}\right&space;),\textup{Beta}\left&space;(&space;y_{2};\alpha_{2},\sum_{&space;k=3}^{m}\alpha&space;_{k}\right&space;),...,\textup{Beta}\left&space;(&space;y_{m-1};\alpha_{m-1},\alpha_{m}\right&space;)$

and the Dirichlet outcome $\mathbf{x}$ is found to be

$\mathbf{x}=\left&space;[&space;x_{1},x_{2},&space;..,x_{k},...x_{m-1},r_{m-1}&space;\right&space;]$

Figure A1 - Left: the q-q plot of $\textup{Beta}\left&space;(&space;x_{1}=y_{1};2,4&space;\right&space;)$. Right: 200 samples of 3D Dirichlet distribution $\textup{Dir}\left&space;(&space;x_{1},x_{2},x_{3};2,2,2&space;\right&space;)$ lying on the plane $x_{1}+x_{2}+x_{3}=1$.

###### Appendix: nucleic acid sequence

A categorical sequence in biology is the nucleic acid sequence (either DNA or RNA), where either acid is made by a chain of linked nucleotides, each one made by a phosphate group, a sugar and one of four canonical/primary nucleo(nitrogenous)bases, adenine (A), cytosin (C), guanine (G), thymine (T), which differentiate the four nucleotides. Thymine (T) is replaced by Uracil (U) in RNA. A sequence of three nucleotides (in total $4^{3}=64$ combinations), known as codon, is the genetic code for synthesizing a protein (long chains of amino acids). Protein coding frames are opened by a start codon, usually the triple AUG, and end with a stop codon. In turn, codons are the core of the protein-coding genes, complex DNA portions capable of activating and regulating the production of the encoded protein. Mutations of genes may result in genetic disorder and relevant diseases.

Alignment of amino acid sequence between a target 3D protein structure and experimental homologous structures is an essential tool for extrapolating and checking the 3D structure of proteins from amino acid sequence. Protein structures are conservative, but not one-to-one related to detected sequences due to variable 3D organization of backbone hydrogen bonds.

The categorical RV is $X=\left&space;\{&space;x_{1}&space;=A,&space;x_{2}=C,x_{3}=G,&space;x_{4}=T/U,x_{5}=Z\right&space;\}$, with m=5. where Z denotes a gap which may be artificially inserted to improve alignment and to keep a uniform length. A typical sequence analysis is the multiple sequence alignment (MSA), in which r sequences of equal length n are aligned to discover conservation or mutation of the nucleotides at a sequence location i=1,…,n. The r aligned sequences form an $r\times&space;n$ matrix, where no column made of gaps (Z) exists. Columns of high purity, i.e. dominated by a single nucleobasis (denoted by a letter), are cornerstones for finding the best alignment. The interest is the distribution of the four nucleotides in a generic column (location) i=1,…,n. The occurrence probability of the category $x_{k}$ in a generic column i is a binomial PMF. To this end, the RV $X_{i}$ becomes a random vector $\mathbf{X}_{i}=\left&space;\{&space;X_{ji}&space;,j=1,...,r\right&space;\}$. Let us denote the number of occurrences (RV) of $x_{k}$ with $N_{k}\left&space;(&space;i&space;\right&space;)=\sum_{j=1}^{r}\left&space;[&space;X_{ji}=x_{k}&space;\right&space;]$. We write the PMF

$P\left&space;(&space;N_{k(i)}=n_{k}(i)&space;\right&space;)=\begin{pmatrix}r\\&space;n_{k}(i)\end{pmatrix}p_{k}(i)^{n_{k}(i)}(1-p_{k}(i))^{r-n_{k}(i)}$

which becomes the likelihood function of an observed sequence $x_{ji}$. The MLE holds $\hat{p}_{k}(i)=n_{k}(i)/r$. The argument $\hat{k}(i)=\textup{arg}\textup{max}_{k}\hat{p}_{k}(i)$ denotes the dominant category (nucleobasis) of the column i. The sequence $\hat{K}=\left&space;\{&space;\hat{k}&space;(i)\right&space;\}$ is the sequence of the column MLE, known as the consensus sequence. Let us denote the MLE of $\hat{k}(i)=\textup{arg}\textup{max}_{k}\hat{p}_{k}(i)$ with $p_{\hat{k}(i)}$, namely $p_{\hat{k}(i)}=\textup{max}_{k}\hat{p}_{k}(i)$. Given a suitable threshold $p_{min}<1$ but close to 1, a consensus sequence such that $p_{\hat{k}(i)}>p_{min}$ is a candidate to declare the alignment acceptable. The frequency of the four nucleotides in protein-coding human genes is fairly uniform, $p_{k}\cong&space;1/4$. It becomes non uniform in non-coding genes.

Figure A1 shows 16 sequences of 16 columns and their column counts. The dominant counts define the consensus sequence CAAATCCGGTTCGAG.

Figure A2 - Left, the counts of each category (nucleobasis) and the max column count. Right: 16 category sequences. The last column is pure.