## 2 EM Algorithm for Maximum Likelihood Reconstruction in Emission Tomography  
The expectation-maximization (EM) algorithm [27] is a statistical algorithm for iteratively computing maximum likelihood estimates when data contain random missing values. Here \"random\" means these missing values do not provide extra information about the parameters we wish to estimate. We first give a brief summary of the EM algorithm below. Since there exist the missing and the observed (or incomplete) components, we can define the complete data set as a combination of the incomplete and the missing data. Note, however, that our aim is to estimate the unknown parameters by maximizing the log-likelihood of the incomplete data. The rationale for the EM algorithm is that if maximizing the incomplete data likelihood is difficult while maximizing the complete data likelihood is easy, then EM can be used to compute iteratively the maximum of the incomplete data likelihood by maximizing the complete data likelihood in each iteration. Let \(\mathcal{C}\) be the complete data set given by \(\mathcal{C}=[\mathcal{Y},\mathcal{M}]\), where \(\mathcal{Y}\) denotes the incomplete data and \(\mathcal{M}\) the missing data. Let \(l_{\mathcal{C}}(\mathbf{x})\) be the log-likelihood based on the complete data \(\mathcal{C}\) and \(l(\mathbf{x})\) the log-likelihood of the incomplete data \(\mathcal{Y}\), where \(\mathbf{x}\) is a \(p\)-vector for the unknown parameters. Let \(\hat{\mathbf{x}}\) be the maximum likelihood (ML) estimate of \(\mathbf{x}\). Then iteration \(k+1\) of the EM algorithm comprises two steps:  
1. **E-Step**: Compute the conditional expectation of the complete data log-likelihood given the incomplete data and \(\mathbf{x}^{(k)}\), and denote this function by \[Q(\mathbf{x};\mathbf{x}^{(k)})=E(l_{\mathcal{C}}(\mathbf{x})\mid\mathcal{Y},\mathbf{x}^{(k)})\] (11)
2. **M-Step**: Update the \(\mathbf{x}\) estimate by maximizing the \(Q\) function, namely \[\mathbf{x}^{(k+1)}=\operatorname*{argmax}_{\mathbf{x}}Q(\mathbf{x};\mathbf{x}^{(k)})\] (12)  
One major advantage of EM is that it guarantees, under certain regularity conditions, that the incomplete data log-likelihood \(l(\mathbf{x})\) increases in consecutive iterations before convergence. Note that EM requires availability of the \(Q\) function in a closed form; otherwise, a Monte-Carlo E-step can be used to replace the E-step [28]. The EM algorithm was first applied to emission tomograph by [PERSON] and [PERSON] [8] and [PERSON] and [PERSON] [29]. Both papers adopt the Poisson model for emission counts, namely \(y_{i}\) are independent Poisson random variables with mean \(\mu_{i}=A_{i}\mathbf{x}\). This model assumes \(r_{i}=0\); otherwise, we can depict \(y_{i}\) as the value after subtracting \(r_{i}\) from the bin \(i\) measurement. From this Poisson model, we can formulate the complete data as \(\mathcal{C}=\{y_{ij}:\;y_{i}=\sum_{j=1}^{p}y_{ij}\}\), where \(y_{ij}\) follows the Poisson distribution with mean \(\mu_{ij}=a_{ij}x_{j}\). Clearly, each \(y_{ij}\) represents the unknown portion of measurement on camera bin \(i\) attributed to image pixel \(j\). The corresponding complete data log-likelihood is  
\[l_{\mathcal{C}}(\mathbf{x})=\sum_{i=1}^{n}\sum_{j=1}^{p}\left\{-\mu_{ij}+y_{ij} \log\mu_{ij}\right\} \tag{13}\]  
and the corresponding \(Q\) function is  
\[Q(\mathbf{x};\mathbf{x}^{(k)})=\sum_{i=1}^{n}\sum_{j=1}^{p}\left\{-\mu_{ij}+y_{ij}^{( k)}\log\mu_{ij}\right\} \tag{14}\]  
where \(y_{ij}^{(k)}=E(y_{ij}\mid y_{i},\mathbf{x}^{(k)})\). Since the conditional distribution of \(y_{ij}\mid y_{i}\) is \(\text{Binomial}(y_{i};\mu_{ij}/\mu_{i})\), we have \(y_{ij}^{(k)}=y_{i}a_{ij}x_{j}^{(k)}/\sum_{t=1}^{p}a_{it}x_{t}^{(k)}\). Thus after solving \(\
abla_{j}Q(\mathbf{x};\mathbf{x}^{(k)})=0\), the M-step of the EM algorithm gives the following updating formula for \(\mathbf{x}\):  
\[x_{j}^{(k+1)}=\frac{x_{j}^{(k)}}{\sum_{i=1}^{n}a_{ij}}\sum_{i=1}^{n}\frac{a_{ ij}y_{i}}{\sum_{t=1}^{p}a_{it}x_{t}^{(k)}} \tag{15}\]  
for \(j=1,\ldots,p\).