Expectation maximization (EM) algorithm is a popular and powerful mathematical method for parameter estimation in case that there exist both observed data and hidden data. The EM process depends on an implicit relationship between observed data and hidden data which is specified by a mapping function in traditional EM and a joint probability density function (PDF) in practical EM. However, the mapping function is vague and impractical whereas the joint PDF is not easy to be defined because of heterogeneity between observed data and hidden data. The research aims to improve competency of EM by making it more feasible and easier to be specified, which removes the vagueness. Therefore, the research proposes an assumption that observed data is the combination of hidden data which is realized as an analytic function where data points are numerical. In other words, observed points are supposedly calculated from hidden points via regression model. Mathematical computations and proofs indicate feasibility and clearness of the proposed method which can be considered as an extension of EM.
Method Article
Expectation Maximization Algorithm with Combinatorial Assumption
https://doi.org/10.21203/rs.3.rs-1304827/v1
This work is licensed under a CC BY 4.0 License
Version 1
posted
You are reading this latest preprint version
expectation maximization (EM)
observed data
hidden data
mapping function
joint probability density function
combinatorial assumption
regression model
Expectation maximization (EM) algorithm developed by Dempster, Laird, and Rubin called DLR (Dempster, Laird, & Rubin, 1977) is an extension of maximum likelihood estimation (MLE) method when there exist both observed data X and hidden data Y and there is an implicit mapping φ: X → Y such that φ–1(Y) = {
: φ(X) = Y}. Let f(X | Θ) be the density probability function (PDF) of random variable X and let g(Y | Θ) be the PDF of random variable Y. Note, Θ is (vector) parameter.
$$g\left(Y|{\Theta }\right)=\underset{{\phi }^{-1}\left(Y\right)}{\int }f\left(X|{\Theta }\right)\text{d}X$$
The conditional PDF of X given Y, denoted k(X | Y, Θ), is specified as follows:
$$k\left(X|Y,{\Theta }\right)=\frac{f\left(X|{\Theta }\right)}{g\left(Y|{\Theta }\right)}$$
Given sample \(\mathcal{Y}\) = {Y1, Y2,…, YN} whose all Yi (s) are mutually independent and identically distributed (iid), EM has many iterations and each iteration has two steps in which expectation step called E-step determines the conditional expectation Q(Θ | Θ(t)) and maximization step (M-step) re-estimates parameter as follows:
E-step:
The expectation Q(Θ | Θ(t)) is determined based on current parameter Θ(t) (Nguyen, 2020).
\(Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)=\sum _{i=1}^{N}\underset{{\phi }^{-1}\left({Y}_{i}\right)}{\int }k\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\text{l}\text{o}\text{g}\left(f\left(X|{\Theta }\right)\right)\text{d}X\) |
(1.1) |
M-step:
The next parameter Θ(t+1) is a maximizer of Q(Θ | Θ(t)) with subject to Θ. Note that Θ(t+1) will become current parameter at the next iteration (the (t+1)th iteration).
EM converges at some tth iteration; at that time, Θ* = Θ(t+1) = Θ(t) is the optimal estimate of EM process. The expectation Q(Θ | Θ(t)) based on the mapping φ(X) represents the traditional EM proposed by DLR. Because it is too vague to specify the mapping φ(X), practical EM issues the joint PDF of X and Y denoted f(X, Y | Θ) as prerequisite condition to run EM such that:
$$f\left(Y|{\Theta }\right)=\underset{X}{\int }f\left(X,Y|{\Theta }\right)\text{d}X$$
$$f\left(X|Y,{\Theta }\right)=\frac{f\left(X,Y|{\Theta }\right)}{\underset{X}{\int }f\left(X,Y|{\Theta }\right)\text{d}X}$$
The expectation Q(Θ | Θ(t)) becomes (Nguyen, 2020):
|
\(Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)=\sum _{i=1}^{N}\underset{X}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\text{l}\text{o}\text{g}\left(f\left(X,{Y}_{i}|{\Theta }\right)\right)\text{d}X\) |
(1.2) |
However, if X and Y are too different in context due to data heterogeneity and they are not independent, it will be difficult to define the PDF f(X, Y | Θ). In general, it is not easy to specify both the mapping φ(X) in traditional EM and the joint PDF f(X, Y | Θ) in practical EM. Therefore, the purpose of this research is to extend competency of EM, in which the observed data Y is assumed to be the combination of X, called combinatorial assumption (CA). In other words, Y can be calculated by an analytic function of X, which is more feasible than specifying the mapping φ(X) and easier than specifying the joint PDF f(X, Y). The analytic function is arbitrary but it is linear called regression function in this research for convenience and feasibility. In some cases, it is possible to convert some analytic functions into linear functions. The next section is the main one which focuses on EM with CA. In general, this research is an extension of EM algorithm.
Almost there is no work to modify EM algorithm itself by support of regression model as methodological extension of EM but many researches utilize EM to learn regression model. Although the context of my methodological research is different, it is necessary to survey some other works related to both EM algorithm and regression model. Zang et al. (Zhang, Deng, & Su, 2016) used EM to learn linear regression model from partial missing data. Their experimental results showed EM algorithm is better than other estimation methods in their case. The essence of regression model estimation is to estimate error random variable with suppose that the error variable distributes normally. If there is no assumption of normal error distribution, EM algorithm may not be effective, which is a drawback of EM. Therefore, Barazandeh et al. (Barazandeh, Ghafelebashi, Razaviyayn, & Sriharsha, 2021) proposed a new algorithm based on combining the alternating direction method of multipliers (ADMM) with EM ideology in order to overcome this drawback. Kokic (Kokic, 2002) also used EM algorithm to learn multivariate regression model in case of missing data and applied the method into time series analysis. Zhang and Rockette (Zhang & Rockette, 2006) proposed an EM algorithm for finding the SML estimate and for variance estimation. Essentially, they utilized EM in case of missing data. Because Gaussian mixture is often learned by EM algorithm, the mixture model of regression functions is also learned by EM algorithm. Faria and Soromenho (Faria & Soromenho, 2009) made a comparison of three methods such as tradition EM, classification EM, and stochastic EM in learning such mixture of regression functions. Their experimental results showed that which method is better is dependent on the true regression line and the initialization step of EM. Griffith (Griffith & Smith, 2010)researched some simplifications of EM on regression model, especially, in case of missing data. In general, most of other researches related to EM and regression models focus on two problems: 1) Learning regression model in case of missing data by EM algorithm because EM is appropriate to impute missing values as hidden data. 2) Learning mixture model of regression functions because mixture model is a popular application of EM.
Given sample \(\mathcal{Y}\) = {y1, y2,…, yN} of size N in which all Yi (s) are mutually independent and identically distributed (iid). Suppose X = (x1, x2,…, xn) which is vector of size n distributes normally with mean vector µand covariance matrix Σ as follows:
|
\(f\left(X|\mu ,{\Sigma }\right)={\left(2\pi \right)}^{-\frac{n}{2}}{\left|{\Sigma }\right|}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}{\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}(X-\mu )\right)\) |
(2.1) |
Where the superscript “T” denotes vector (matrix) transposition operator. The nx1 mean µ and the nxn covariance matrix Σ:
$$\mu =\left(\begin{array}{c}{\mu }_{1}\\ {\mu }_{2}\\ ⋮\\ {\mu }_{n}\end{array}\right),{\Sigma }=\left(\begin{array}{cccc}{\sigma }_{1}^{2}& {\sigma }_{12}^{2}& \cdots & {\sigma }_{1n}^{2}\\ {\sigma }_{21}^{2}& {\sigma }_{2}^{2}& \cdots & {\sigma }_{2n}^{2}\\ ⋮& ⋮& \ddots & ⋮\\ {\sigma }_{n1}^{2}& {\sigma }_{n2}^{2}& \cdots & {\sigma }_{n}^{2}\end{array}\right)$$
Note, σij2 is covariance of xi and xj whereas σi2 is variance of xi. Let Y = (y1, y2,…, ym)T be random variable representing all sample random variable Yi = (yi1, yi2,…, yim)T. Note, X is nx1 vector and Y is mx1 vector. Suppose there is an assumption that Y is a combination of partial random variables (components) of X such that:
$${y}_{i}={\alpha }_{i0}+\sum _{j=1}^{n}{\alpha }_{ij}{x}_{i}$$
As a generalization, let A be mxn matrix whose elements are called regressive coefficients as follows:
\(A=\left(\begin{array}{ccccc}{\alpha }_{10}& {\alpha }_{11}& {\alpha }_{12}& \cdots & {\alpha }_{1n}\\ {\alpha }_{20}& {\alpha }_{21}& {\alpha }_{22}& \cdots & {\alpha }_{2n}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{m0}& {\alpha }_{m1}& {\alpha }_{m2}& \cdots & {\alpha }_{mn}\end{array}\right)\) \({A}_{0}=\left(\begin{array}{c}{\alpha }_{10}\\ {\alpha }_{20}\\ ⋮\\ {\alpha }_{m0}\end{array}\right),\tilde{A}=\left(\begin{array}{cccc}{\alpha }_{11}& {\alpha }_{12}& \cdots & {\alpha }_{1n}\\ {\alpha }_{21}& {\alpha }_{22}& \cdots & {\alpha }_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {\alpha }_{m1}& {\alpha }_{m2}& \cdots & {\alpha }_{mn}\end{array}\right)\)
This implies
|
\(Y={A}_{0}+\tilde{A}X\) |
(2.2) |
As a convention, let:
$${\tilde{A}}_{i}={\left({\alpha }_{i1},{\alpha }_{i2},\dots ,{\alpha }_{in}\right)}^{T}$$
Then
$$\tilde{A}=\left(\begin{array}{c}{\tilde{A}}_{1}\\ {\tilde{A}}_{2}\\ ⋮\\ {\tilde{A}}_{m}\end{array}\right)$$
$${y}_{i}={\alpha }_{i0}+{\tilde{A}}_{i}^{T}X$$
The equation above is regression function in which Y is called responsor and X is called regressor whereas A is called regressive matrix. The assumption is combinatorial assumption (CA) aforementioned and the method proposed here is called CA method or CA algorithm. Suppose Y distributes normally with mean \({A}_{0}+\tilde{A}X\) and covariance matrix Sas follows:
|
\(f\left(Y|X,A,S\right)={\left(2\pi \right)}^{-\frac{m}{2}}{\left|S\right|}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}{\left(Y-{A}_{0}-\tilde{A}X\right)}^{T}{S}^{-1}(Y-{A}_{0}-\tilde{A}X)\right)\) |
(2.3) |
The mxm covariance matrix S is:
$$S=\left(\begin{array}{cccc}{s}_{1}^{2}& {s}_{12}^{2}& \cdots & {s}_{1m}^{2}\\ {s}_{21}^{2}& {s}_{2}^{2}& \cdots & {s}_{2m}^{2}\\ ⋮& ⋮& \ddots & ⋮\\ {s}_{m1}^{2}& {s}_{m2}^{2}& \cdots & {s}_{m}^{2}\end{array}\right)$$
Note, sij2 is covariance of yi and yj whereas si2 is variance of yi. As a convention, let:
$${S}_{i}={\left({s}_{i1},{s}_{i2},\dots ,{s}_{in}\right)}^{T}$$
Then
$$S=\left(\begin{array}{c}{S}_{1}^{T}\\ {S}_{2}^{T}\\ ⋮\\ {S}_{n}^{T}\end{array}\right)=\left({S}_{1},{S}_{2},\dots ,{S}_{n}\right)$$
As a convention, \({\stackrel{-}{S}}_{i}\) is the one which is like Si except that S is replaced by S−1. Note, \({\stackrel{-}{S}}_{i}\) is vector. Similarly, notations \({\stackrel{-}{s}}_{ij}^{2}\) and \({\stackrel{-}{s}}_{i}^{2}\) imply such a meaning. The marginal PDF of Y is now defined by support of regression model as follows:
$$f\left(Y|{\Theta }\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(X,Y|{\Theta }\right)\text{d}X\stackrel{\scriptscriptstyle\text{def}}{=}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|\mu ,{\Sigma }\right)f\left(Y|X,A,S\right)\text{d}X$$
Where parameter Θ = (µ, Σ, A, S)T is compound parameter. The equation above is not real total probability rule but it implies that the conditional PDF f(Y | X) is substituted by regression model. Consequently, the expectation Q(Θ | Θ(t)) becomes:
|
\(Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\text{l}\text{o}\text{g}\left(f\left(X|\mu ,{\Sigma }\right)f\left({Y}_{i}|X,A,S\right)\right)\text{d}X\) |
(2.4) |
It is necessary to specify the conditional PDF f(X | Y, Θ). Indeed, we have:
$$f\left(X|Y,{\Theta }\right)=\frac{f\left(X|\mu ,{\Sigma }\right)f\left(Y|X,A,S\right)}{\underset{{\mathbb{R}}^{n}}{\int }f\left(X|\mu ,{\Sigma }\right)f\left(Y|X,A,S\right)\text{d}X}=\frac{f\left(X,y|{\Theta }\right)}{f\left(y|{\Theta }\right)}$$
The joint PDF f(X, Y | Θ) which is the numerator of f(X | Y, Θ) is defined:
$$f\left(X,Y|{\Theta }\right)\stackrel{\scriptscriptstyle\text{def}}{=}f\left(X|\mu ,{\Sigma }\right)f\left(Y|X,A,S\right)={\left(2\pi \right)}^{-\frac{n+m}{2}}{\left(\left|{\Sigma }\right|\left|S\right|\right)}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\left(Y-{A}_{0}\right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\mu \right)\right)*{f}_{0}\left(X,Y|{\Theta }\right)$$
Where,
$${f}_{0}\left(X,y|{\Theta }\right)=\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)+{\left(\tilde{A}\mu \right)}^{T}{S}^{-1}\tilde{A}\mu \right)\right)$$
The expression \({\left(\tilde{A}X\right)}^{T}{S}^{-1}\tilde{A}X\) is approximated with µ as follows:
$${\left(\tilde{A}X\right)}^{T}{S}^{-1}\tilde{A}X\cong {\left(\tilde{A}\mu \right)}^{T}{S}^{-1}\tilde{A}\mu$$
As a result, f0(X, Y | Θ) and f(X, Y | Θ) is approximated as follows:
$${f}_{0}\left(X,Y|{\Theta }\right)\cong \text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\right)\right)$$
And
$$f\left(X,Y|{\Theta }\right)\cong {\left(2\pi \right)}^{-\frac{n+m}{2}}{\left(\left|{\Sigma }\right|\left|S\right|\right)}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(Y-{A}_{0}-\tilde{A}\mu \right)}^{T}{S}^{-1}\left(Y-{A}_{0}-\tilde{A}\mu \right)\right)\right)*{f}_{0}\left(X,y|{\Theta }\right)$$
The approximation by removing X-dependency from the expression \({\left(\tilde{A}X\right)}^{T}{S}^{-1}\tilde{A}X\) is reasonable because the PDF f(X | Y, Θ) contains second-order proportion with the built-in expression \({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)\) and this PDF also reflects regression model with another built-in expression \({\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\) including parameter S2. In other words, the dependency of \({\left(\tilde{A}X\right)}^{T}{S}^{-1}\tilde{A}X\) on X is unnecessary. Moreover, EM process will adjust parameters by the best way later. In following proofs and computations, we will see that such dependency removal also makes the research easy to apply shifted Gaussian integral.
The denominator of f(X | Y, Θ) which is f(Y | Θ) is the integral of f(X, Y | Θ) over X:
$$f\left(Y|{\Theta }\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(X,Y|{\Theta }\right)\text{d}X={\left(2\pi \right)}^{-\frac{n+m}{2}}{\left(\left|{\Sigma }\right|\left|S\right|\right)}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(Y-{A}_{0}-\tilde{A}\mu \right)}^{T}{S}^{-1}\left(Y-{A}_{0}-\tilde{A}\mu \right)\right)\right)B$$
Where B is the integral of f0(X, Y | Θ) over X:
$$B=\underset{{\mathbb{R}}^{n}}{\int }{f}_{0}\left(X,Y|{\Theta }\right)\text{d}X\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\right)\right)\text{d}X$$
It requires to calculate B to determine f(X | Y, Θ). By referring the appendix, we can denote:
$$B\cong E\left({f}_{0}\left(X,y|{\Theta }\right)\right)$$
Obviously, B is totally determined. Thus, f(Y | Θ) is approximated as follows:
$$f\left(Y|{\Theta }\right)\cong {\left(2\pi \right)}^{-\frac{n+m}{2}}{\left(\left|{\Sigma }\right|\left|S\right|\right)}^{-\frac{1}{2}}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(Y-{A}_{0}-\tilde{A}\mu \right)}^{T}{S}^{-1}\left(Y-{A}_{0}-\tilde{A}\mu \right)\right)\right)*E\left({f}_{0}\left(X,y|{\Theta }\right)\right)$$
As a result, the PDF f(X | Y, Θ) is approximated as follows:
$$f\left(X|Y,{\Theta }\right)=\frac{f\left(X,Y|{\Theta }\right)}{f\left(Y|{\Theta }\right)}\cong \frac{1}{E\left({f}_{0}\left(X,Y|{\Theta }\right)\right)}\text{*}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\right)\right)$$
Let k(Y|Θ) be the constant with subject to X but it is function of Ywith parameter Θ, which is defined as follows:
|
\(k\left(Y|{\Theta }\right)=\frac{1}{E\left({f}_{0}\left(X,Y|{\Theta }\right)\right)}\) |
(2.5) |
Shortly, the conditional PDF f(X | Y, Θ(t)) is specified (approximated) at E-step of some tthiteration process as follows:
|
\(f\left(X|Y,{{\Theta }}^{\left(t\right)}\right)\stackrel{\scriptscriptstyle\text{def}}{=}k\left(Y|{{\Theta }}^{\left(t\right)}\right)\text{*}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\right)\right)\) |
(2.6) |
Consequently, the expectation Q(Θ | Θ(t)) at E-step of some tth iteration is totally determined. At M-step of current tth iteration, Q(Θ|Θ(t)) is maximized by setting its partial derivatives regarding Θ to be zero. The first-order partial derivative of Q(Θ | Θ(t)) with regard to µ with note that Q(Θ | Θ(t)) is analytic function is:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \mu }=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\frac{\partial \text{l}\text{o}\text{g}\left(f\left(X|\mu ,{\Sigma }\right)\right)}{\partial \mu }\text{d}X$$
$$=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\text{d}X$$
$$=\left(\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\left(X-{\mu }^{\left(t\right)}\right)}^{T}\text{d}X-N{\mu }^{T}+N{\left({\mu }^{\left(t\right)}\right)}^{T}\right){{\Sigma }}^{-1}$$
By referring to the appendix, we have:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \mu }=\left(-N{\mu }^{T}+N{\left({\mu }^{\left(t\right)}\right)}^{T}+\sum _{i=1}^{N}E\left(X|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)\right){{\Sigma }}^{-1}$$
Note, Σ is invertible and symmetric. As a convention, I denote:
$$E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)=E\left(X|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)$$
The next parameter µ(t+1) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is solution of the equation \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \mu }={0}^{T}\), with note that 0 = (0, 0,…, 0)T is zero vector, as follows:
|
\({\mu }^{\left(t+1\right)}={\mu }^{\left(t\right)}+\frac{1}{N}\sum _{i=1}^{N}E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\) |
(2.7) |
The first-order partial derivative of Q(Θ | Θ(t)) with regard to Σ with note that Q(Θ | Θ(t)) is analytic function is:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {\Sigma }}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\frac{\partial \text{l}\text{o}\text{g}\left(f\left(X|\mu ,{\Sigma }\right)\right)}{\partial {\Sigma }}\text{d}X$$
$$=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(-\frac{1}{2}{{\Sigma }}^{-1}+\frac{1}{2}{{\Sigma }}^{-1}\left(X-\mu \right){\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\right)\text{d}X$$
The next parameter Σ(t+1) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is solution of equation formed by setting \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {\Sigma }}\) to zero matrix. Let (0) denote zero matrix.
$$\left(0\right)=\left(\begin{array}{cccc}0& 0& \cdots & 0\\ 0& 0& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & 0\end{array}\right)$$
We have:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {\Sigma }}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(-{\Sigma }+\left(X-{\mu }^{\left(t\right)}\right){\left(X-{\mu }^{\left(t\right)}\right)}^{T}\right)\text{d}X=\left(0\right)$$
Note, µ is replaced by µ(t). Thus, the next parameter Σ(t+1) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is obtained:
$${{\Sigma }}^{\left(t+1\right)}=\frac{1}{N}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(X-{\mu }^{\left(t\right)}\right){\left(X-{\mu }^{\left(t\right)}\right)}^{T}\text{d}X$$
By referring to the appendix, we obtain:
|
\({{\Sigma }}^{\left(t+1\right)}=\frac{1}{N}\sum _{i=1}^{N}E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\) |
(2.8) |
As a convention, I denote:
$$E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)=E\left(X{X}^{T}|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)$$
The first-order partial derivative of Q(Θ | Θ(t)) with regard to A0 with note that Q(Θ | Θ(t)) is analytic function is:
$${\left(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {A}_{0}}\right)}^{T}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\frac{\partial \text{l}\text{o}\text{g}\left(f\left(X|Y,{\Theta }\right)\right)}{\partial {\alpha }_{0}}\text{d}X$$
$$={S}^{-1}\left(\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({Y}_{i}-{A}_{0}-\tilde{A}{\mu }^{\left(t\right)}-\tilde{A}\left(X-{\mu }^{\left(t\right)}\right)\right)\text{d}X\right)$$
$$={S}^{-1}\left(\sum _{i=1}^{N}{Y}_{i}-N{A}_{0}-N\tilde{A}{\mu }^{\left(t\right)}-\tilde{A}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(X-{\mu }^{\left(t\right)}\right)\text{d}X\right)$$
By referring to the appendix, we obtain:
$${\left(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {A}_{0}}\right)}^{T}={S}^{-1}\left(\sum _{i=1}^{N}{Y}_{i}-N{A}_{0}-N{\tilde{A}}^{\left(t\right)}{\mu }^{\left(t\right)}-{\tilde{A}}^{\left(t\right)}\sum _{i=1}^{N}E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)$$
Note, \(\tilde{A}\) is replaced by \({\tilde{A}}^{\left(t\right)}\) at current tth iteration. Therefore, the next parameter A0(t+1) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is obtained by setting the partial derivative \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial {A}_{0}}\) to be zero:
|
\({A}_{0}^{\left(t+1\right)}=\frac{1}{N}\left(\sum _{i=1}^{N}{Y}_{i}-N{\tilde{A}}^{\left(t\right)}{\mu }^{\left(t\right)}-N{\tilde{A}}^{\left(t\right)}\left({\mu }^{\left(t+1\right)}-{\mu }^{\left(t\right)}\right)\right)\) |
(2.9) |
Note, the next parameter µ(t+1) is specified by equation 2.7. The first-order partial derivative of Q(Θ | Θ(t)) with regard to \(\tilde{\alpha }\) with note that Q(Θ | Θ(t)) is analytic function is (Saliba, 2016):
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \tilde{A}}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\frac{\partial \text{l}\text{o}\text{g}\left(f\left(X|{Y}_{i},{\Theta }\right)\right)}{\partial \tilde{A}}\text{d}X$$
$$=2\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\left({Y}_{i}-{A}_{0}-\tilde{A}X\right)}^{T}{S}^{-1}\left({X}^{T}⨂{I}_{m}\right)\text{d}X$$
Where \(⨂\) denotes Kronecker product and Im is mxm identity matrix. The next parameter \({\tilde{A}}^{\left(t+1\right)}\) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is obtained by setting the partial derivative \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \tilde{A}}\) to be zero. Because the Kronecker product \({X}^{T}⨂{I}_{n}\) occurs, the equation \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial \tilde{A}}\) being equal to zero is equivalent to 1 equation \({u}_{j\ne 1}\left(\tilde{A}\right)=0\) and n equations \({v}_{j}\left(\tilde{A}\right)=0\), where
$${u}_{j\ne 1}\left(\tilde{A}\right)={\stackrel{-}{S}}_{j}^{T}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({Y}_{i}-{A}_{0}-\tilde{A}X\right)\text{d}X$$
$$={\stackrel{-}{S}}_{j}^{T}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({Y}_{i}-{A}_{0}-\tilde{A}\left(X-\mu \right)-\tilde{A}\mu \right)\text{d}X$$
$$={\stackrel{-}{S}}_{j}^{T}\left(\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)-N\tilde{A}\left(\mu +E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)\right)$$
And
$${v}_{j}\left(\tilde{A}\right)={\stackrel{-}{S}}_{1}^{T}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({Y}_{i}-{A}_{0}-\tilde{A}X\right){x}_{j}\text{d}X$$
$$={\stackrel{-}{S}}_{1}^{T}\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(\left({Y}_{i}-{A}_{0}\right)\left({x}_{j}-{\mu }_{j}\right)+\left({Y}_{i}-{A}_{0}\right){\mu }_{j}-\tilde{A}X\right){x}_{j}\text{d}X$$
$$={\stackrel{-}{S}}_{1}^{T}\left(\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\left(E\left({x}_{j}|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)+{\mu }_{j}\right)-\tilde{A}E\left({x}_{j}X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)$$
$$={\stackrel{-}{S}}_{1}^{T}\left(\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\left(E\left({x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)+{\mu }_{j}\right)-N\tilde{A}\stackrel{-}{E}\left({x}_{j}X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)$$
Where,
$$\stackrel{-}{E}\left({x}_{j}X|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)=\left(\begin{array}{c}\stackrel{-}{E}\left({x}_{1}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\\ \stackrel{-}{E}\left({x}_{2}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\\ ⋮\\ \stackrel{-}{E}\left({x}_{n}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\end{array}\right)$$
It is necessary to determine every partial expectation \(\stackrel{-}{E}\left({x}_{i}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\) to determine the expectation \(\stackrel{-}{E}\left({x}_{j}X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\). Indeed, by referring to the appendix, we have:
$$\stackrel{-}{E}\left({x}_{i}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)=\stackrel{-}{E}\left(\left({x}_{i}-{\mu }_{i}^{\left(t\right)}\right)\left({x}_{j}-{\mu }_{j}^{\left(t\right)}\right)+{x}_{i}{\mu }_{j}+{x}_{j}{\mu }_{i}^{\left(t\right)}-{\mu }_{i}^{\left(t\right)}{\mu }_{j}^{\left(t\right)}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)$$
$$=E\left({x}_{i}{x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)+{\mu }_{i}^{\left(t\right)}E\left({x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)+{\mu }_{j}^{\left(t\right)}E\left({x}_{i}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)-{\mu }_{i}^{\left(t\right)}{\mu }_{j}^{\left(t\right)}$$
As a result, the next parameter \({\tilde{A}}^{\left(t+1\right)}\) is solution of n+1 equations as follows:
|
\({\tilde{A}}^{\left(t+1\right)}\stackrel{\scriptscriptstyle\text{def}}{=}\left\{\begin{array}{l}\tilde{A}\left({\mu }^{\left(t\right)}+E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)=\frac{1}{N}\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\\ {\left.\tilde{A}\stackrel{-}{E}\left({x}_{j}X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)=\frac{1}{N}\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\left(E\left({x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)+{\mu }_{j}^{\left(t\right)}\right)\right|}_{j=\stackrel{-}{1,n}}\end{array}\right.\) |
(2.10) |
Note, µ is replaced by the current µ(t). The equation 2.10 can be solved by Newton-Raphson method. The first-order partial derivative of Q(Θ | Θ(t)) with regard to S with note that Q(Θ | Θ(t)) is analytic function is:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial S}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\frac{\partial \text{l}\text{o}\text{g}\left(f\left(X|{Y}_{i},{\Theta }\right)\right)}{\partial S}\text{d}X$$
$$=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(-\frac{1}{2}{S}^{-1}+\frac{1}{2}{S}^{-1}\left({Y}_{i}-{A}_{0}-\tilde{A}X\right){\left({Y}_{i}-{A}_{0}-\tilde{A}X\right)}^{T}{S}^{-1}\right)\text{d}X$$
The next parameter S(t+1) at M-step of some tth iteration that maximizes Q(Θ|Θ(t)) is solution of equation formed by setting \(\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial S}\) to zero matrix. By referring to the appendix, we have:
$$\frac{\partial Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)}{\partial S}=\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(-S+\left({Y}_{i}-{A}_{0}-\tilde{A}X\right){\left({Y}_{i}-{A}_{0}-\tilde{A}X\right)}^{T}\right)\text{d}X$$
$$=-NS+\sum _{i=1}^{N}\underset{{\mathbb{R}}^{n}}{\int }f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left(\left({Y}_{i}-{A}_{0}\right){\left({Y}_{i}-{A}_{0}\right)}^{T}-2\tilde{A}\left(X-\mu \right){\left({Y}_{i}-{A}_{0}\right)}^{T}-2\tilde{A}\mu {\left({Y}_{i}-{A}_{0}\right)}^{T}+\tilde{A}\left(X-\mu \right){\left(X-\mu \right)}^{T}{\tilde{A}}^{T}+2\tilde{A}\left(X-\mu \right){\mu }^{T}{\tilde{A}}^{T}+\tilde{A}\mu {\mu }^{T}{\tilde{A}}^{T}\right)\text{d}X$$
$$=-NS+\tilde{A}\mu {\mu }^{T}{\tilde{A}}^{T}+\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}-2\tilde{A}\mu \right){\left({Y}_{i}-{A}_{0}\right)}^{T}+\sum _{i=1}^{N}\left(-2\tilde{A}E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({\left({Y}_{i}-{A}_{0}\right)}^{T}-{\mu }^{T}{\tilde{A}}^{T}\right)+\tilde{A}E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\tilde{A}}^{T}\right)=\left(0\right)$$
Therefore, we obtain:
|
\({S}^{\left(t+1\right)}=\frac{1}{N}\left({\tilde{A}}^{\left(t\right)}\mu {\mu }^{T}{\left({\tilde{A}}^{\left(t\right)}\right)}^{T}+\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}^{\left(t\right)}-2{\tilde{A}}^{\left(t\right)}\mu \right){\left({Y}_{i}-{A}_{0}^{\left(t\right)}\right)}^{T}+\sum _{i=1}^{N}\left(-2{\tilde{A}}^{\left(t\right)}E\left(X|f\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)\left({\left({Y}_{i}-{A}_{0}^{\left(t\right)}\right)}^{T}-{\mu }^{T}{\left({\tilde{A}}^{\left(t\right)}\right)}^{T}\right)+{\tilde{A}}^{\left(t\right)}E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\left({\tilde{A}}^{\left(t\right)}\right)}^{T}\right)\right)\) |
(2.11) |
Note, µ, A0 and \(\tilde{A}\) are replaced by µ(t), A0(t) and \({\tilde{A}}^{\left(t\right)}\), respectively. In general, CA method is EM process with two steps as follows:
E-step:
Determining the conditional PDF f(X | Y, Θ(t)) specified by equation 2.6 based on current parameter Θ(t).
$$f\left(X|Y,{{\Theta }}^{\left(t\right)}\right)\stackrel{\scriptscriptstyle\text{def}}{=}k\left(Y|{{\Theta }}^{\left(t\right)}\right)\text{*}\text{e}\text{x}\text{p}\left(-\frac{1}{2}\left({\left(X-\mu \right)}^{T}{{\Sigma }}^{-1}\left(X-\mu \right)-2{\left(Y-{A}_{0}\right)}^{T}{S}^{-1}\tilde{A}\left(X-\mu \right)\right)\right)$$
M-step:
Calculating next parameters Θ(t+1) = (µ(t+1), Σ(t+1), A(t+1), S(t+1))T based on f(X | Y, Θ(t)) determined in the E-step, specified by equations 2.7, 2.8, 2.9, 2.10, and 2.11.
|
\({\mu }^{\left(t+1\right)}={\mu }^{\left(t\right)}+\frac{1}{N}\sum _{i=1}^{N}E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\) \({{\Sigma }}^{\left(t+1\right)}=\frac{1}{N}\sum _{i=1}^{N}E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\) \({A}_{0}^{\left(t+1\right)}=\frac{1}{N}\left(\sum _{i=1}^{N}{Y}_{i}-N{\tilde{A}}^{\left(t\right)}{\mu }^{\left(t\right)}-N{\tilde{A}}^{\left(t\right)}\left({\mu }^{\left(t+1\right)}-{\mu }^{\left(t\right)}\right)\right)\) \({\tilde{A}}^{\left(t+1\right)}\stackrel{\scriptscriptstyle\text{def}}{=}\left\{\begin{array}{l}\tilde{A}\left(\mu +E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\right)=\frac{1}{N}\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\\ {\left.\tilde{A}\stackrel{-}{E}\left({x}_{j}X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)=\frac{1}{N}\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}\right)\left(E\left({x}_{j}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)+{\mu }_{j}\right)\right|}_{j=\stackrel{-}{1,n}}\end{array}\right.\) \({S}^{\left(t+1\right)}=\frac{1}{N}\left({\tilde{A}}^{\left(t\right)}\mu {\mu }^{T}{\left({\tilde{A}}^{\left(t\right)}\right)}^{T}+\sum _{i=1}^{N}\left({Y}_{i}-{A}_{0}^{\left(t\right)}-2{\tilde{A}}^{\left(t\right)}\mu \right){\left({Y}_{i}-{A}_{0}^{\left(t\right)}\right)}^{T}+\sum _{i=1}^{N}\left(-2{\tilde{A}}^{\left(t\right)}E\left(X|{Y}_{i},{{\Theta }}^{\left(t\right)}\right)\left({\left({Y}_{i}-{A}_{0}^{\left(t\right)}\right)}^{T}-{\mu }^{T}{\left({\tilde{A}}^{\left(t\right)}\right)}^{T}\right)+{\tilde{A}}^{\left(t\right)}E\left(X{X}^{T}|{Y}_{i},{{\Theta }}^{\left(t\right)}\right){\left({\tilde{A}}^{\left(t\right)}\right)}^{T}\right)\right)\) |
In practice, it is not necessary to compute the covariance matrix Σ(t+1) and the variance S(t+1) because computational cost is high and it is also really ineffective to estimate Σ and S. Note, the condition that both Σ and S are invertible along with |Σ| ≠ 0 and |S| ≠ 0 is not easy to assert over many computational iterations. The most important parameters are µ and A and we should fix the other parameters Σ and S with hints of predefined bias or background knowledge.
Although combinatorial assumption (CA) is subjective assumption, it will reach high usefulness and high effectiveness if there is strong regressive relationship between observed data and hidden data in many cases. The regression function may not be linear but it is easy to convert nonlinear functions into linear function in some cases. For instance, there are some nonlinear functions as follows:
| Logarithm function | \(y={\alpha }_{0}+\sum _{j=1}^{n}{\alpha }_{j}\text{l}\text{o}\text{g}\left({x}_{j}\right)\) |
|---|---|
| Exponent function | \(y=\text{e}\text{x}\text{p}\left({\alpha }_{0}+\sum _{j=1}^{n}{\alpha }_{j}\text{l}\text{o}\text{g}{x}_{j}\right)\) |
| Product function | \(y={\alpha }_{0}\prod _{j=1}^{n}{\left({x}_{j}\right)}^{{\alpha }_{j}}\) |
Given logarithm function, the transformation is \(y={\alpha }_{0}+{\tilde{\alpha }}^{T}U\) and uj = log(xj). For exponent function, the transformation is \(v={\alpha }_{0}+{\tilde{\alpha }}^{T}X\) and v = log(y). For product function, the transformation is \(v={\beta }_{0}+{\tilde{\alpha }}^{T}U\), v = log(y), uj = log(xj), and β0 = log(α0).
Recall that the dimensions of X and Y are n and m, respectively. Note, if m ≥ n, there is no information loss which is ideal case of CA method. Otherwise, if m < n, the important parameters such as µ and A vary much with large amplitude due to information loss of dimension reduction. The problem is called large-scale variation. Therefore, in practice there should be restrictions on µ and A specified as constraints u(µ) ≤ 0 and v(A) ≤ 0; for example, slope of regressive hyperplane specified by the normal vectors \(-{\tilde{A}}_{i}\), which varies in a predefined interval, is a constraint. The solution is to apply Lagrange duality method into maximizing Q(Θ|Θ(t)), in which a Lagrange function l(Θ, κ, λ | Θ(t)) is defined as sum of Q(Θ|Θ(t)) and constraints u(µ) ≤ 0 and v(A) ≤ 0 as follows:
$$l\left({\Theta },\kappa ,{\lambda }|{{\Theta }}^{\left(t\right)}\right)=Q\left({\Theta }|{{\Theta }}^{\left(t\right)}\right)+\kappa u\left(\mu \right)+{\lambda }v\left(A\right)$$
Note Q(Θ|Θ(t)) is maximized indirectly by maximizing the Lagrange function l(Θ, κ, λ | Θ(t)), in which κ ≥ 0 and λ ≥ 0 called Lagrange multipliers are concerned. Another simple trick to alleviate the large-scale variation of µ and A is to initialize appropriate values µ(0) and A(0) at the first iteration of EM process.
An application of CA is to learn Bayesian parameter. According to Bayesian statistics, the parameter Θ is considered as random variable. Given sample D = {X1, X2,…, XN) whose all observations Xi (s) are iid, the posterior probability of Θ given D denoted P(Θ | D, ξ) is calculated based on D and the prior probability of Θ denoted P(Θ | ξ).
$$P\left({\Theta }|D,\xi \right)=\frac{P\left(D|{\Theta },\xi \right)P\left({\Theta }|\xi \right)}{P\left(D|\xi \right)}$$
Where ξ denotes background knowledge about Θ which can be considered sub-parameter but we do not focus on learning ξ. Let X be random variable representing all Xi. The most popular method to learn Bayesian parameter is to use binomial sample and beta distribution. Here CA method considers the parameter Θ as hidden data and random variable X as observed data.
$$X={\alpha }_{0}+{\tilde{\alpha }}^{T}{\Theta }$$
Therefore, the posterior probability P(Θ | D, ξ) is the same to the conditional PDF f(X | y, Θ) aforementioned. After EM process is finished, f(X | y, Θ) or P(Θ | D, ξ) is obviously determined.
In general, CA is an extension of EM, whose strong point is feasibility and simplicity with clearness of mathematic formulas. Especially, in the ideal case that dimension of Y (y) is larger than or equal to dimension of X, CA method will result out best estimates because there is no information loss. However, its drawback is the large-scale variation in determining important parameters µ and A when dimension of Y (y) is smaller than dimension of X. Therefore, in the future I will research Lagrange duality method to maximize the expectation Q(Θ|Θ(t)) with constraints on µ and A.
This appendix focuses on calculating integrals related to the following Gaussian function:
$$f\left(x\right)=\text{e}\text{x}\text{p}\left(-{\left(Ax\right)}^{T}R\left(Ax\right)-{\left(Bx\right)}^{T}S\left(Bx\right)+{c}^{T}x\right)$$
Note, x = (x1, x2,…, xh)T is hx1 vector, A is nxh matrix, B is mxh matrix, and c is hx1 vector whereas R is nxn matrix and S is mxm matrix. Moreover, R and S are invertible and positive definite. Firstly, we calculate the integral of f(x).
$$\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)\text{d}x=\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{\left(Ax\right)}^{T}R\left(Ax\right)-{\left(Bx\right)}^{T}S\left(Bx\right)+{c}^{T}x\right)dx$$
Diagonalizing R and S as follows:
$$R=U{\Gamma }{U}^{T}$$
$$V=V{\Lambda }{V}^{T}$$
Where nxn matrix U and mxm matrix V are orthonormal and nxn matrix Γ and mxm matrix Λ are eigenmatrices:
$$U=\left(\begin{array}{cccc}{u}_{11}& {u}_{12}& \cdots & {u}_{1n}\\ {u}_{21}& {u}_{22}& \cdots & {u}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {u}_{n1}& {u}_{n2}& \cdots & {u}_{nn}\end{array}\right),V=\left(\begin{array}{cccc}{v}_{11}& {v}_{12}& \cdots & {v}_{1n}\\ {v}_{21}& {v}_{22}& \cdots & {v}_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ {v}_{n1}& {v}_{n2}& \cdots & {v}_{nn}\end{array}\right)$$
$${\Gamma }=\left(\begin{array}{cccc}{\gamma }_{1}& 0& \cdots & 0\\ 0& {\gamma }_{2}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\gamma }_{n}\end{array}\right),{\Lambda }=\left(\begin{array}{cccc}{\lambda }_{1}& 0& \cdots & 0\\ 0& {\lambda }_{2}& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\lambda }_{n}\end{array}\right)$$
We have (Mathprof, 2013):
$$\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)\text{d}x=\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{\left(Ax\right)}^{T}U{\Gamma }{U}^{T}\left(Ax\right)-{\left(Bx\right)}^{T}V{\Lambda }{V}^{T}\left(Bx\right)+{c}^{T}x\right)dx$$
Let,
$${y}_{i}={\gamma }_{i}\sum _{j=1}^{h}{x}_{j}\sum _{k=1}^{n}{a}_{kj}{u}_{ki},{z}_{i}={\lambda }_{i}\sum _{j=1}^{h}{x}_{j}\sum _{k=1}^{n}{b}_{kj}{v}_{ki}$$
We have:
$$\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)\text{d}x=\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-\sum _{i=1}^{n}{y}_{t}^{2}-\sum _{i=1}^{m}{z}_{t}^{2}+\sum _{i=1}^{h}{c}_{t}{x}_{t}\right)dx$$
Let
$${\stackrel{-}{a}}_{ij}=\sum _{k=1}^{n}{a}_{kj}{u}_{ki},{\stackrel{-}{b}}_{ij}=\sum _{i=1}^{m}{b}_{kj}{v}_{ki}$$
This implies
$${y}_{i}^{2}={\gamma }_{i}^{2}{\left(\sum _{j=1}^{h}{\stackrel{-}{a}}_{ij}{x}_{j}\right)}^{2}={\gamma }_{t}^{2}\left(\sum _{j=1}^{h}{\stackrel{-}{a}}_{ij}^{2}{x}_{j}^{2}+2\sum _{j=1}^{h}\sum _{k=1,k\ne j}^{h}{\stackrel{-}{a}}_{ij}{\stackrel{-}{a}}_{ik}{x}_{j}{x}_{k}+\sum _{j=1}^{h}{\stackrel{-}{a}}_{ij}^{2}\right)$$
$${z}_{i}^{2}={\lambda }_{i}^{2}{\left(\sum _{j=1}^{h}{\stackrel{-}{b}}_{ij}{x}_{j}\right)}^{2}={\lambda }_{i}^{2}\left(\sum _{j=1}^{h}{\stackrel{-}{b}}_{ij}^{2}{x}_{j}^{2}+2\sum _{j=1}^{h}\sum _{k=1,k\ne j}^{h}{\stackrel{-}{b}}_{ij}{\stackrel{-}{b}}_{ik}{x}_{j}{x}_{k}+\sum _{j=1}^{h}{\stackrel{-}{b}}_{ij}^{2}\right)$$
Suppose X distributes normally with mean µ = (µ1, µ2,…, µh)T, the product xjxl is approximated by µjµl and then, the quantities yi2 and zi2 are re-written as follows:
$${y}_{i}^{2}\cong {\gamma }_{i}^{2}\left(\sum _{j=1}^{h}{\stackrel{-}{a}}_{ij}^{2}{x}_{j}^{2}+\sum _{j=1}^{h}{\mu }_{j}\sum _{k=1,k\ne j}^{h}\left(2{\stackrel{-}{a}}_{ij}{\stackrel{-}{a}}_{ik}{\mu }_{k}+{\stackrel{-}{a}}_{ij}^{2}\right)\right)$$
$${z}_{i}^{2}\cong {\lambda }_{i}^{2}\left(\sum _{j=1}^{h}{\stackrel{-}{b}}_{ij}^{2}{x}_{j}^{2}+\sum _{j=1}^{h}{\mu }_{j}\sum _{k=1,k\ne j}^{h}\left(2{\stackrel{-}{b}}_{ij}{\stackrel{-}{b}}_{ik}{\mu }_{k}+{\stackrel{-}{b}}_{ij}^{2}\right)\right)$$
Let
$${p}_{ij}={\mu }_{j}\sum _{k=1,k\ne j}^{h}\left(2{\stackrel{-}{a}}_{ij}{\stackrel{-}{a}}_{ik}{\mu }_{k}+{\stackrel{-}{a}}_{ij}^{2}\right)$$
$${q}_{ij}={\mu }_{j}\sum _{k=1,k\ne j}^{h}\left(2{\stackrel{-}{b}}_{ij}{\stackrel{-}{b}}_{ik}{\mu }_{k}+{\stackrel{-}{b}}_{ij}^{2}\right)$$
We have:
$${y}_{i}^{2}\cong {\gamma }_{i}^{2}\sum _{j=1}^{h}\left({\stackrel{-}{a}}_{ij}^{2}{x}_{j}^{2}+{p}_{ij}\right)$$
$${z}_{i}^{2}\cong {\lambda }_{i}^{2}\sum _{j=1}^{h}\left({\stackrel{-}{b}}_{ij}^{2}{x}_{j}^{2}+{q}_{ij}\right)$$
The integral \(\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)dx\) is re-written as follows:
$$\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)dx\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-\sum _{i=1}^{n}\sum _{j=1}^{h}{\gamma }_{i}^{2}\left({\stackrel{-}{a}}_{ij}^{2}{x}_{j}^{2}+{p}_{ij}\right)-\sum _{i=1}^{m}\sum _{j=1}^{h}{\lambda }_{i}^{2}\left({\stackrel{-}{b}}_{ij}^{2}{x}_{j}^{2}+{q}_{ij}\right)+\sum _{j=1}^{h}{x}_{j}{c}_{j}\right)dx$$
Let w be the maximum value among n, m, and h and let
As a convention, let
$$E\left(f\left(x|{\Theta }\right)\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right)dx$$
$$E\left({x}_{i}|f\left(x|{\Theta }\right)\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right){x}_{i}dx$$
$$E\left({x}_{i}^{2}|f\left(x|{\Theta }\right)\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right){x}_{i}^{2}dx$$
$$E\left({x}_{i}{x}_{j}|f\left(x|{\Theta }\right)\right)=\underset{{\mathbb{R}}^{n}}{\int }f\left(x\right){x}_{i}{x}_{j}dx$$
Where Θ implies parameters of f(x) such as A, R, B, S, c. We have (Wikipedia, 2021):
$$E\left(f\left(x|{\Theta }\right)\right)\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(\sum _{i=1}^{w}\left(-{F}_{i}{x}_{i}^{2}+{G}_{i}{x}_{i}+{H}_{i}\right)\right)dx$$
$$=\prod _{i=1}^{w}\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{F}_{i}{x}_{i}^{2}+{G}_{i}{x}_{i}+{H}_{i}\right)dx=\prod _{i=1}^{w}{\left(\frac{\pi }{{F}_{i}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)$$
Note, if Fi = 0, the respective expression \({\left(\frac{\pi }{{F}_{i}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)\) is removed from the product series. We calculate the integral E(xi | f(x|Θ)) as follows (Wikipedia, 2021):
$$E\left({x}_{i}|f\left(x|{\Theta }\right)\right)\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(\sum _{j=1}^{w}\left(-{F}_{j}{x}_{j}^{2}+{G}_{j}{x}_{j}+{H}_{j}\right)\right){x}_{i}dx$$
$$=\left(\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{F}_{i}{x}_{i}^{2}+{G}_{i}{x}_{i}+{H}_{i}\right){x}_{i}d{x}_{i}\right)*\prod _{j=1,j\ne i}^{w}{\left(\frac{\pi }{{F}_{j}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{j}^{2}}{4{F}_{j}}+{H}_{j}\right)$$
$$=\frac{\sqrt{\pi }{G}_{i}}{2{F}_{i}^{3/2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)*\prod _{j=1,j\ne i}^{w}{\left(\frac{\pi }{{F}_{j}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{j}^{2}}{4{F}_{j}}+{H}_{j}\right)$$
Consequently, the integral E(xi2 | f(x|Θ)) as follows (Wikipedia, 2021):
$$E\left({x}_{i}^{2}|f\left(x|{\Theta }\right)\right)\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(\sum _{j=1}^{w}\left(-{F}_{j}{x}_{j}^{2}+{G}_{j}{x}_{j}+{H}_{j}\right)\right){x}_{i}^{2}dx$$
$$=\left(\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{F}_{i}{x}_{i}^{2}+{G}_{i}{x}_{i}+{H}_{i}\right){x}_{i}^{2}d{x}_{i}\right)*\prod _{j=1,j\ne i}^{w}{\left(\frac{\pi }{{F}_{j}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{j}^{2}}{4{F}_{j}}+{H}_{j}\right)$$
$$=\frac{\sqrt{\pi }\left(2{F}_{i}+{G}_{i}^{2}\right)}{4{F}_{i}^{5/2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)*\prod _{j=1,j\ne i}^{w}{\left(\frac{\pi }{{F}_{j}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{j}^{2}}{4{F}_{j}}+{H}_{j}\right)$$
Moreover, the integral E(xixj | f(x|Θ)) as follows:
$$E\left({x}_{i}{x}_{j}|f\left(x|{\Theta }\right)\right)\cong \underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(\sum _{k=1}^{w}\left(-{F}_{k}{x}_{k}^{2}+{G}_{k}{x}_{k}+{H}_{k}\right)\right){x}_{i}{x}_{j}dx$$
$$=\left(\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{F}_{i}{x}_{i}^{2}+{G}_{i}{x}_{i}+{H}_{i}\right){x}_{i}d{x}_{i}\right)*\left(\underset{{\mathbb{R}}^{n}}{\int }\text{e}\text{x}\text{p}\left(-{F}_{j}{x}_{j}^{2}+{G}_{j}{x}_{j}+{H}_{j}\right){x}_{j}d{x}_{j}\right)*\prod _{k=1,k\ne i,k\ne j}^{w}{\left(\frac{\pi }{{F}_{k}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{k}^{2}}{4{F}_{k}}+{H}_{k}\right)$$
$$=E\left({x}_{i}|f\left(x|{\Theta }\right)\right)E\left({x}_{j}|f\left(x|{\Theta }\right)\right)*\prod _{k=1,k\ne i,k\ne j}^{w}{\left(\frac{\pi }{{F}_{k}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{k}^{2}}{4{F}_{k}}+{H}_{k}\right)$$
In general, we have:
\(E\left(f\left(x|{\Theta }\right)\right)\cong \prod _{i=1}^{w}{\left(\frac{\pi }{{F}_{i}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)\) \(E\left({x}_{i}|f\left(x|{\Theta }\right)\right)\cong \frac{\sqrt{\pi }\left(2{F}_{i}+{G}_{i}^{2}\right)}{4{F}_{i}^{5/2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{i}^{2}}{4{F}_{i}}+{H}_{i}\right)*\prod _{j=1,j\ne i}^{w}{\left(\frac{\pi }{{F}_{j}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{j}^{2}}{4{F}_{j}}+{H}_{j}\right)\) \(E\left({x}_{i}{x}_{j}|f\left(x|{\Theta }\right)\right)\cong E\left({x}_{i}|f\left(x|{\Theta }\right)\right)E\left({x}_{j}|f\left(x|{\Theta }\right)\right)*\prod _{k=1,k\ne i,k\ne j}^{w}{\left(\frac{\pi }{{F}_{k}}\right)}^{\frac{1}{2}}\text{e}\text{x}\text{p}\left(\frac{{G}_{k}^{2}}{4{F}_{k}}+{H}_{k}\right)\)
By extending these integrals with vector x, we determine:
\(E\left(x|f\left(x|{\Theta }\right)\right)=\left(\begin{array}{c}E\left({x}_{1}|f\left(x|{\Theta }\right)\right)\\ E\left({x}_{2}|f\left(x|{\Theta }\right)\right)\\ ⋮\\ E\left({x}_{h}|f\left(x|{\Theta }\right)\right)\end{array}\right)\) \(E\left(x{x}^{T}|f\left(x|{\Theta }\right)\right)=\left(\begin{array}{cccc}E\left({x}_{1}^{2}|f\left(x|{\Theta }\right)\right)& E\left({x}_{1}{x}_{2}|f\left(x|{\Theta }\right)\right)& \cdots & E\left({x}_{1}{x}_{h}|f\left(x|{\Theta }\right)\right)\\ E\left({x}_{2}{x}_{1}|f\left(x|{\Theta }\right)\right)& E\left({x}_{2}^{2}|f\left(x|{\Theta }\right)\right)& \cdots & E\left({x}_{2}{x}_{h}|f\left(x|{\Theta }\right)\right)\\ ⋮& ⋮& \ddots & ⋮\\ E\left({x}_{h}{x}_{1}|f\left(x|{\Theta }\right)\right)& E\left({x}_{h}{x}_{2}|f\left(x|{\Theta }\right)\right)& \cdots & E\left({x}_{h}^{2}|f\left(x|{\Theta }\right)\right)\end{array}\right)\)
Barazandeh, B., Ghafelebashi, A., Razaviyayn, M., & Sriharsha, R. (2021, May 12). Efficient Algorithms for Estimating the Parameters of Mixed Linear Regression Models. arXiv. Retrieved from https://arxiv.org/pdf/2105.05953
Dempster, A. P., Laird, N. M., & Rubin, D. B. (1977). Maximum Likelihood from Incomplete Data via the EM Algorithm. (M. Stone, Ed.) Journal of the Royal Statistical Society, Series B (Methodological), 39(1), 1-38.
Faria, S., & Soromenho, G. (2009, February 11). Fitting mixtures of linear regressions. Journal of Statistical Computation and Simulation, 80(2), 201-225. doi:10.1080/00949650802590261
Griffith, D. A., & Smith, A. (2010). Some simplifications for the expectationmaximization (em) algorithm: the linear regression model case. InterStat, 23. doi:10.1.1.613.8033
Kokic, P. (2002). The EM Algorithm for a Multivariate Regression Model: including its applications to a non-parametric regression model and a multivariate time series model. University of York, Department of Computer Science. York: University of York. Retrieved from https://www.cs.york.ac.uk/euredit/_temp/The%20Euredit%20Software/NAG%20Prototype%20platform/WorkingPaper4.pdf
Mathprof. (2013, March 22). Multidimensional Gaussian integral. (PlanetMath ) Retrieved November 23, 2021, from PlanetMath: https://planetmath.org/multidimensionalgaussianintegral
Nguyen, L. (2020). Tutorial on EM algorithm. MDPI. Preprints. doi:10.20944/preprints201802.0131.v8
Saliba, G. (2016, July 21). Not understanding derivative of a matrix-matrix product. (Stack Overflow) Retrieved December 06, 2021, from Mathematics Stack Exchange: https://math.stackexchange.com/questions/1866757/not-understanding-derivative-of-a-matrix-matrix-product
Wikipedia. (2021, December 3). List of integrals of exponential functions. (Wikimedia Foundation) Retrieved December 6, 2021, from Wikipedia website: https://en.wikipedia.org/wiki/List_of_integrals_of_exponential_functions
Zhang, X., Deng, J., & Su, R. (2016). The EM algorithm for a linear regression model with application to a diabetes data. 2016 International Conference on Progress in Informatics and Computing (PIC). Shanghai: IEEE. doi:10.1109/PIC.2016.7949477
Zhang, Z., & Rockette, H. E. (2006, November 4). An EM algorithm for regression analysis with incomplete covariate information. Journal of Statistical Computation and Simulation. doi:10.1080/10629360600565202
Competing interests:
The authors declare no competing interests.