04-numanal.Rmd

# Numerical Methods

## Introduction

### Computer representation of real numbers

Any positive decimal number $x$ is represented by the ordered coefficents $\{ d_j : j = n, n - 1, \ldots \} \subseteq \{ 0, 1, \ldots, 9 \}$

\begin{equation}
  x = d_n 10^n + d_{n - 1} 10 ^{n - 1} + \cdots + d_1 10 + d_0 + d_{-1} 10^{-1} + \cdots
  (\#eq:decimal)
\end{equation}

For same number $x$, other base $2$ can also be used with binary digits $\{ a_j \} \subseteq \{ 0, 1 \}$

\begin{equation}
  x = a_k 2^k + a_{k - 1} 2 ^{k - 1} + \cdots + a_1 2 + a_0 + a_{-1} 2^{-1} + \cdots
  (\#eq:base2)
\end{equation}

Point between $a_0$ and $a_{-1}$ is called the radix point here.

```{r}
sfsmisc::digitsBase(320, base = 10)
sfsmisc::digitsBase(320, base = 2)
```

See Equations \@ref(eq:decimal) and \@ref(eq:base2). Numbers are expressed with series.

```{example, machineeps, name = "Identical and nearly equal"}
$0.3 - 0.1$ is equal to $0.2$. Can we check this?
```

```{r}
(.3 - .1) == .2
```

It is obviously same, but `R` says it is different. Why?

```{r}
.Machine$double.eps
```

The above number is the smallest positive floating number that the machine can recognize. `all.equal()` function can solve this kind of near-equality problem.

```{r}
all.equal(.2, .3 - .1)
```


## Root-finding in One Dimension

In statistics, it is one of issues to find solutions of

$$f(x) = 0$$

There are various algorithms.

### Bisection method

```{r bisecfig, echo=FALSE, fig.cap="Illustration of bisection method"}
tibble(x = c(0, 2)) %>% 
  ggplot(aes(x = x)) +
  modelr::geom_ref_line(h = 0) +
  stat_function(fun = function(x) {
    (x - 3)^2 - 2
  }) +
  geom_vline(xintercept = c(.5, 1.7), alpha = .5, col = I("grey10")) +
  annotate(
    geom = "point",
    x = -sqrt(2) + 3,
    y = 0,
    col = "red",
    shape = 1
  )
```

Figure \@ref(fig:bisecfig) presents the motivation of bisection method. On both sides of the root, one side of function value is positive and the other side is negative. Thus, if we find any set like this, then we only narrow the two points until finding the solution.

\begin{algorithm}[H] \label{alg:bisection}
  \SetAlgoLined
  \SetKwInOut{Input}{input}
  \SetKwInOut{Output}{output}
  \Input{Equation system $f(x) = 0$, error bound $\epsilon$}
  Initialize two points $x_0$ and $x_1$ such that $$f(x_0) f(x_1) \le 0$$\;
  \If{$f(x_0) f(x_1) < 0$}{
    Change initial values\;
  }
  Set error $e = \lvert x_1 - x_0 \rvert$\;
  \While{$e > \epsilon$}{ \label{alg:bisecwhile}
    Half $$x_2 = \frac{x_0 + x_1}{2}$$\;
    Length of the interval becomes half $e = \frac{e}{2}$\;
    \eIf{$f(x_0) f(x_2) < 0$}{
      Update $x_1 = x_2$\;
    }{
      Update $x_0 = x_2$\;
    }
  }
  \Output{$x = x_2$}
  \caption{Bisection algorithm}
\end{algorithm}

In Line $\ref{alg:bisecwhile}$, we can use condition

$$\lvert f(x_2) \rvert > \epsilon$$

instead, which means that we did not find the root yet.

```{example, bisecexm}
Solve $$a^2 + y^2 + \frac{2ay}{n - 1} = n - 2$$

where $a$ is a specified constant and $n > 2$ is an integer.
```

```{solution}
It can be shown that the analytical solution is

$$y = - \frac{a}{n - 1} \pm \sqrt{n - 2 + a^2 + \Big( \frac{a}{n - 1} \Big)^2}$$
```

```{r}
f_bisec <- function(x, a = .5, n = 20) {
  a^2 + x^2 + 2 * a * x / (n - 1) - (n - 2)
}
#-----------------------------
bisection <- function(x0, x1, fun, eps = .Machine$double.eps^.25, rep_max = 1000, ...) {
  iter <- 0 # stop too many iteration
  if (fun(x0, ...) * fun(x1, ...) > 0) {
    stop(gettextf("both %s and %s should be satisfy the condition", expression(x0), expression(x1)))
  }
  init <- seq(x0, x1, length.out = 3) # x0 x2 x1
  y <- f_bisec(init)
  while (iter < 1000 && abs(y[2]) > eps) {
    iter <- iter + 1
    if (y[1] * y[2] < 0) {
      init[3] <- init[2]
      y[3] <- y[2]
    } else {
      init[1] <- init[2]
      y[1] <- y[2]
    }
    init[2] <- (init[1] + init[3]) / 2
    y[2] <- fun(init[2], ...)
  }
  c(init[2], y[2])
}
```

Using initioal values $x_0 = 0$ and $X_1 = 100$,

```{r}
(bi_exm <- bisection(0, 100, fun = f_bisec, a = .5, n = 20))
```

$x = `r bi_exm[1]`$ has been computed. The following figure shows that this answer is reasonable.

```{r bisecroot, echo=FALSE, fig.cap="Example curve"}
tibble(x = c(0, 10)) %>% 
  ggplot(aes(x = x)) +
  modelr::geom_ref_line(h = 0) +
  stat_function(fun = f_bisec) +
  annotate(
    geom = "point", 
    x = bi_exm[1], 
    y = f_bisec(bi_exm[1]), 
    col = "red", 
    shape = 4,
    size = 3
  )
```

### Brent's method

Brent's method combines the root bracketing and bisection with inverse quadratic interpolation. `uniroot()` uses this method. Refer to Example \@ref(exm:bisecexm).

```{r}
(brent <- 
  uniroot(
    f = f_bisec,
    interval = c(0, 100),
    a = .5,
    n = 20
  ))
```

This method assures convergence of the bisection method. Morover, it is generally faster than bisection.


## Numerical Integration

Try to compute

$$I = \int_a^b f(x) dx$$

### Trapezoidal rule

From definition of Riemann integration, we can compute integration $I$ by partitioning intervals. Areas of rectangles can be considered or trapezoids can also be considered. If we use trapezoids, it will be more closed to the target curve, but the formula might be quite complex. For the length of subintervals $h = \frac{b - a}{n}$,

\begin{equation}
  \frac{h}{2} f(a) + h \sum_{i = 1}^{n - 1} f(x_i) + \frac{h}{2} f(b)
  (\#eq:trape)
\end{equation}

For fun, we use `Rcpp` for trapezoid method. `Rcpp` integrate `R` and `C++`. This accelerate execution speed like loop.

```{r, eval=FALSE}
library(Rcpp)
```

The following code should be written in `cpp` file separately, or in `cppFunction()` as character.

```{Rcpp}
#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector trapezoid(Function target, double a, double b, int n) {
  double h = (b - a) / n;
  NumericVector fa = target(a);
  NumericVector fb = target(b);
  NumericVector integral = (fa + fb) / 2;
  double x = a;
  NumericVector fx = target(x);
  
  for(int i = 0; i < n; i++) {
    x += h;
    NumericVector fx = target(x);
    integral += fx;
  }
  
  integral = integral * h;
  
  return(integral);
}
```

Consider standard normal densitiy. Compare

$$P(-1.96 \le Z \le 1.96)$$

```{r}
phi <- function(x) {
  1 / sqrt(2 * pi) * exp(- x^2 / 2)
}
#----------------------------------
tibble(x1 = -1.96, x2 = 1.96) %>% 
  summarise(
    trapezoid = trapezoid(
      phi,
      a = x1,
      b = x2,
      n = 100
    ),
    pnorm = pnorm(x2) - pnorm(x1)
  )
```

### Adaptive quadrature method

`R` provides a function `integrate()`. This implements a method called an *adaptive quadrature method*. Get

$$\int_0^{\infty} \frac{1}{(\cosh y - \rho r)^{n - 1}}dy$$

with $\rho \in (-1, 1)$, $r \in (-1, 1)$, and $n \ge 2$ integer.

```{r}
integrate_exm <- function(y, n, r, rho) {
  (cosh(y) - rho * r)^(1 - n)
}
```

Denote that $\rho$, $r$, and $n$ should be pre-specified. Consider $(0.2, 0.5, 10)$.

```{r}
integrate(
  f = integrate_exm,
  lower = 0,
  upper = Inf,
  n = 10,
  r = .5,
  rho = .2
)
```


## Maximum Likelihood Problems

Maximum likelihood estimator (MLE) is a estimator such that maximizes likelihood function. Given random sample $x_1, \ldots, x_n \iid f(x_i ; \theta)$, likelihood function can be given by

$$L(\theta) = \prod_{i = 1}^n f(x_i)$$

Then MLE $\hat\theta$ is

\begin{equation}
  \hat\theta = \argmax_{\theta \in \Theta} L(\theta)
  (\#eq:mle)
\end{equation}

Denote that it is equivalent to maximizing log-likilihood $l(\theta) := \ln L(\theta)$.

\begin{equation}
  \hat\theta = \argmax_{\theta \in \Theta} l(\theta)
  (\#eq:mle2)
\end{equation}

Either for $L$ or $l$, we can find the critical point by differentiating in a mathematical point of view.

$$
\begin{cases}
  \frac{d}{d \theta} l(\theta) = 0 \\
  \frac{d^2}{d \theta^2} l(\theta) > 0
\end{cases}
$$

Ignoring the second line, try to find root of first dervative. Finding MLE becomes *root-finding of first derivative function* problem. What we need are

1. Likelihood function or log-likelihood function
2. Its derivative

```{example, findmle, name = "Exponential distribution"}
Let $Y_1, Y_2 \iid Exp(\theta)$, i.e.

$$f(y) = \theta e^{-\theta y}, \quad y > 0, \theta > 0$$

Then the likelihood function is

$$L(\theta) = \theta^2 e^{-\theta(y_1 + y_2)}$$

and log-likelihood

$$l(\theta) = 2 \ln \theta - \theta(y_1 + y_2)$$

Find its MLE $\hat\theta$.
```

```{solution}
Note that for $\theta > 0$,

$$\frac{d}{d\theta}l(\theta) = \frac{2}{\theta} - (y_1 + y_2)$$

Hence, we know that the analytical solution is

$$\hat\theta = \frac{2}{y_1 + y_2}$$
```

```{r}
y <- c(.043, .502)
```

Give input as $l$ and $(y_1, y_2) = (`r y`)$. Here we will use `D()` which enables to output analytical derivative function for `expression`. For example,

```{r}
D(expression(2 * log(theta) - theta * (y1 + y2)), name = "theta")
```

Then we can make the following function.

```{r}
find_mle <- function(l, args, name, interval = c(1, 5), ...) {
  differ <- D(substitute(l), name = name)
  args[[name]] <- 0
  differ_fun <- function(x) {
    args[[name]] <- x
    eval(differ, envir = args, enclos = parent.frame())
  }
  uniroot(
    f = differ_fun,
    interval = interval,
    ...
  )$root
}
#-------------------------------------
find_mle(
  2 * log(theta) - theta * (y1 + y2), 
  args = list(y1 = y[1], y2 = y[2]), 
  name = "theta",
  interval = c(1, 5)
)
```

In `stats4` library, there is a function called `mle()`. We can also use this one.

```{r}
exp_logLik <- function(theta = 1) {
  - length(y) * log(theta) + theta * sum(y) # -l(theta)
}
#-----------------------------------
stats4::mle(exp_logLik) %>% 
  stats4::summary()
```

## One-Dimensional Optimization

In the last section, our custom function finding MLE have tried to find root. On the contrary, `stats4::mle()` optimizes given negative log-likelihood, i.e. find its minimum.

```{example, maxlike, name = "Find maximum of univariate function"}
Maximize the function

$$f(x) = \frac{\ln (1 + \ln x)}{\ln(1 + x)}$$

with respect to $x$.
```

```{r logfrac, fig.cap="Function $f$ in Example"}
log_frac <- function(x) {
  log(1 + log(x)) / log(1 + x)
}
#---------------------------
gg_curve(
  log_frac, from = 2, to = 14, 
  ylab = expression(log(1 + log(x)) / log(1 + x))
)
```

`nlm()` finds minimization with a Newton-type algorithm. `optimize()` performs optimization based on various `method`. To find maximum, we should specify `maximum = TRUE`. It is set to be `FALSE` by default and find the minimum.

```{r}
optimize(
  log_frac,
  lower = 2,
  upper = 8,
  maximum = TRUE
)
```

`maximum` is a point where maximum is occurred and `objective` is a maximum value of the function.


## Two-Dimensional Optimization

```{example, mixgam}
Let $X_1, \ldots, X_3 \indep f(x \mid \lambda_j) \equiv Gamma(r = \frac{1}{2}, \lambda_j)$, i.e.

$$f(x \mid \lambda_j) = \frac{\lambda^{\frac{1}{2}}}{\Gamma(\frac{1}{2})} x^{- \frac{1}{2}} e^{-\lambda x}$$

Then set a mixture $Y \sim f$ s.t. with mixing probaiblity $p_1, p_2, p_3 \;(p_1 + p_2 + p_3 = 1)$.

$$f(y) = \sum_{j = 1}^3 p_j f_j (y \mid \lambda_j)$$

provided that $\lambda_1 + \lambda_2 + \lambda_3 = 1$.
```

`optim()` can be used for multi-parameter optimization. It finds minimum of given function in `fn`. Since we are interested in maximum, we have to *add minus sign in front of the final result to get minimum*.

$$\max f = - \min ( - f)$$

For simplicity, set `param = c(p1, p2, lambda1, lambda2)`. First compute log-likelihood for mixture.

```{r}
mix_ll <- function(param, y) {
  # mixing probability
  prob <- param[1:2]
  prob <- c(prob, 1 - sum(prob))
  # rate of dgamma
  rate <- param[3:4]
  rate <- c(rate, 1 - sum(rate)) # constraint
  dens <-
    sapply(rate, function(b) {
      dgamma(x = y, shape = 1 / 2, rate = 1 / (2 * b))
    }) %*%
    diag(prob) %>% # p_j * f_j
    rowSums()
  -sum(log(dens)) # negative log-likelihood
}
```

```{r}
thet_lam <- c(.6, .25, .15)
```

Let `r thet_lam` be the true $(\lambda_1, \lambda_2, \lambda_3)^T$. Also, $p_1 = p_2 = p_3 = \frac{1}{3}$.

1. Generate random number from this parameter, which play a role of observed sample
2. Set an initial value for each parameter
3. Each step: Get log-likelihood value from this sample using pre-defined `mix_ll()`
4. Find minimum

Following `y` is used for `mix_ll(y = y)`, i.e. sample for likelihood.

```{r optim_init}
init_lam <- 
  sample( # p1 = p2 = p3 = 1/3
    thet_lam, # true parameter
    size = 2000,
    replace = TRUE
  )
y <- rgamma(2000, shape = 1 / 2, rate =  1 / (2 * init_lam)) # mixture gamma with 1/3
```

In `optim()`, initial value should be supplied to `par`.

```{r}
(opt <- optim(par = c(.3, .3, .5, .3), fn = mix_ll, y = y))
```

`$value` is the minimum value of `fn`. `$par` is our interest, parameter vector where the function is minimized. `broom::tidy()` also has a method for `optim` result. This changes it to `tibble` with two columns - `parameter` (parameter name) and `value` (its value).

```{r}
opt_df <- 
  broom::tidy(opt) %>% 
  spread(parameter, value)
#-------------------------
colnames(opt_df) <- c("prob1", "prob2", "lambda1", "lambda2")
opt_df %>% 
  mutate(
    prob3 = 1 - prob1 - prob2,
    lambda3 = 1 - lambda1 - lambda2
  ) %>% 
  select(order(colnames(.)))
```

Compare this MLE to the true $\lambda_j$ `r thet_lam`.


## EM Algorithm

EM algorithm is often applied to find MLE, especially when data are incomplete. This section is mainly following explanation of @Bilmes:1998tg. Let $\mathcal{X}$ be the observed data set. Assume that this is incomplete. Then let $\mathcal{Y}$ be the missing part. Assume that the complete data set exists $\mathcal{Z} = (\mathcal{X}, \mathcal{Y})$. Then the joint density of $\mathcal{Z}$ is

\begin{equation}
  p(\mathbf{z} \mid \boldsymbol\theta) = p(\mathbf{x}, \mathbf{y} \mid \boldsymbol\theta) = p(\mathbf{y} \mid \mathbf{x}, \boldsymbol\theta) p(\mathbf{x} \mid \boldsymbol\theta)
  (\#eq:completejoint)
\end{equation}

MLE problem tries to find the maximum of

\begin{equation}
  L(\boldsymbol\theta \mid \mathcal{X}, \mathcal{Y}) = p(\mathcal{X}, \mathcal{Y} \mid \boldsymbol\theta)
  (\#eq:completelike)
\end{equation}

which is called *complete-data likelihood*. For this kind of likelihood, we use EM algorithm. **E** stands for expectation and **M** for maximization. EM algorithm consists of these two step.

### Expectation step (E-step)

At first, expectation step finds conditional expectation of complete-data likelihood \@ref(eq:completelike) or its log given observed sample and current parameter estimates. Denote that random $\mathcal{Y}$ has not been observed. So we should remove it. Conditinal expectation does this job.

Let $\boldsymbol\theta^{(i - 1)}$ be the current parameter estimates. EM algorithm keeps updating the parameters. At each step, we would have each updated value. What we want is conditional expectation under $\mathcal{Y} \mid \mathcal{X}, \boldsymbol\theta$. Then we should know its conditional density, i.e.

\begin{equation}
  f(\mathbf{y} \mid \mathcal{X}, \boldsymbol\theta)
  (\#eq:missingcond)
\end{equation}

Then E-step calculate the following. Writing the support of $\mathcal{Y}$ by $\mathbb{Y}$,

\begin{equation}
  Q(\boldsymbol\theta, \boldsymbol\theta^{(i - 1)}) := E \bigg[ \ln p(\mathcal{X}, \mathcal{Y} \mid \boldsymbol\theta) \; \Big\vert \; \mathcal{X}, \boldsymbol\theta^{(i - 1)} \bigg] = \int_{\mathbb{Y}} \ln p(\mathcal{X}, \mathbf{y} \mid \boldsymbol\theta) f(\mathbf{y} \mid \mathcal{X}, \boldsymbol\theta^{(i - 1)}) d \mathbf{y}
  (\#eq:estep)
\end{equation}

### Maximization step (M-step)

Maximization step maximizes the conditional expectation $Q$ \@ref(eq:estep) with respect to $\boldsymbol\theta$ given $\boldsymbol\theta^{(i - 1)}$ [@Rizzo:2007aa].

\begin{equation}
  \boldsymbol\theta^{(i)} = \argmax_{\boldsymbol\theta} Q(\boldsymbol\theta, \boldsymbol\theta^{(i - 1)})
  (\#eq:mstep)
\end{equation}

E-step \@ref(eq:estep) and M-step \@ref(eq:mstep) are repeated. It is guaranteed that the log-likelihood is always increasing and the algorithm converges to local maximum.

### EM algorithm for a mixture model

EM algorithm is widely used for mixture model. Recall that

$$p(\mathbf{x} \mid \boldsymbol\beta) = \sum_{m = 1}^M \alpha_m p_m(\mathbf{x} \mid \theta_m)$$

with parameter vector $\boldsymbol\beta = (\alpha_1, \ldots, \alpha_M, \theta_1, \ldots, \theta_M)$. $\alpha_m$ is an mixing probability satisfying that

$$\sum_{m = 1}^M \alpha_m = 1$$

Each $\theta_m$ parameterizes individual density $p_m$. Consider complete sample $\mathcal{Z} = (\mathcal{X}, \mathcal{Y}) = \{ (\mathbf{x}_i, \mathbf{y}_i) \}_1^N$ with unobserved $\mathcal{Y}$. Assume that for each $i = 1, \ldots, N$,

$$\mathbf{x}_i \sim p_k(\mathbf{x} \mid \theta_k) \Rightarrow y_i = k$$

Given $\mathcal{Y}$ values, the log of complete-data likelihood is

\begin{equation}
  \begin{split}
    l(\boldsymbol\beta \mid \mathcal{X}, \mathcal{Y}) & = \ln p(\mathcal{X}, \mathcal{Y} \mid \boldsymbol\beta) \\
    & = \sum_{i = 1}^N \ln p(\mathbf{x}_i \mid \mathbf{y}_i) p(\mathbf{y}_i) \\
    & = \sum_{i = 1}^N \ln \alpha_{\mathbf{y}_i} p_{\mathbf{y}_i}(\mathbf{x}_i \mid \theta_{\mathbf{y}_i})
  \end{split}
  (\#eq:logmix)
\end{equation}

Recall that $\mathcal{Y}$ has not been observed. So let $\mathcal{Y}$ be random. In each step we update (guess) the parameter for above likelihood \@ref(eq:logmix). Write it as

$$\boldsymbol\beta^g = (\alpha_1^g, \ldots, \alpha_M^g, \theta_1^g, \ldots, \theta_M^g)$$

Given $\boldsymbol\beta^g$, compute $p_j(\mathbf{x}_i \mid \theta_j^g)$ for each $i$ and $j$. Bayes rule implies that

\begin{equation}
  \begin{split}
    p(\mathbf{y}_i \mid \mathbf{x}_i, \boldsymbol\beta^g) & = \frac{\alpha_{\mathbf{y}_i}^g p_{\mathbf{y}_i}(\mathbf{x}_i \mid \theta_{\mathbf{y}_i}^g)}{p_j(\mathbf{x}_i \mid \theta_j^g)} \\
    & = \frac{\alpha_{\mathbf{y}_i}^g p_{\mathbf{y}_i}(\mathbf{x}_i \mid \theta_{\mathbf{y}_i}^g)}{\sum\limits_{m = 1}^M \alpha_m^g p_m(\mathbf{x}_i \mid \theta_m^g)}
  \end{split}
  (\#eq:mixbayes)
\end{equation}

Denote that $\alpha_j = P(\text{component}\; j)$ are kind of prior probabilities. Let $\mathbf{y} = (\mathbf{y}_1, \ldots, \mathbf{y}_N)$. Then E-step \@ref(eq:estep) becomes

\begin{equation}
  \begin{split}
    Q(\boldsymbol\beta, \boldsymbol\beta^g) & = \sum_{\mathbf{y} \in \mathbb{Y}} l(\boldsymbol\beta \mid \mathcal{X}, \mathbf{y}) p(\mathbf{y} \mid \mathcal{X}, \boldsymbol\beta^g) \\
    & = \sum_{\mathbf{y} \in \mathbb{Y}} {\color{blue}{\sum_{i = 1}^N \ln \alpha_{\mathbf{y}_i} p_{\mathbf{y}_i}(\mathbf{x}_i \mid \theta_{\mathbf{y}_i})}} \prod_{j = 1}^N p(\mathbf{y}_j \mid \mathbf{x}_j, \boldsymbol\beta^g) \\
    & = \underbrace{\sum_{y_1 = 1}^M \sum_{y_2 = 1}^M \cdots \sum_{y_N = 1}^M}_{\mathbb{Y}} {\color{blue}{\sum_{i = 1}^N \sum_{l = 1}^M \delta_{l, \mathbf{y}_i} \ln \alpha_l p_l (\mathbf{x}_i \mid \theta_l)}} \prod_{j = 1}^N p(\mathbf{y}_j \mid \mathbf{x}_j, \boldsymbol\beta^g) \\
    & = {\color{blue}{\sum_{l = 1}^M \sum_{i = 1}^N \ln \alpha_l p_l (\mathbf{x}_i \mid \theta_l)}} \underbrace{\sum_{y_1 = 1}^M \sum_{y_2 = 1}^M \cdots \sum_{y_N = 1}^M {\color{blue}{\delta_{l, \mathbf{y}_i}}} \prod_{j = 1}^N p(\mathbf{y}_j \mid \mathbf{x}_j, \boldsymbol\beta^g)}_{(\ast)}
  \end{split}
  (\#eq:mixcond)
\end{equation}

To simplify Equation \@ref(eq:mixcond), see $(\ast)$ part.

\begin{equation}
  \begin{split}
    (\ast) & = \bigg( \sum_{y_1 = 1}^M \cdots \sum_{y_{i - 1} = 1}^M \sum_{y_{i + 1} = 1}^M \cdots \sum_{y_N = 1}^M \prod_{j \neq i}^N p(\mathbf{y}_j \mid \mathbf{x}_j, \boldsymbol\beta^g) \bigg) \\
    & = \prod_{j \neq i}^N \bigg( \underbrace{\sum_{y_j = 1}^M p(\mathbf{y}_j \mid \mathbf{x}_j, \boldsymbol\beta^g)}_{= 1} \bigg) p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) \\
    & = p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) \quad \because \sum_{i = 1}^N p(i \mid \mathbf{x}_j, \boldsymbol\beta^g) = 1
  \end{split}
  (\#eq:mixsimple)
\end{equation}

From Equations \@ref(eq:mixcond) and \@ref(eq:mixsimple), we can conclude E-step.

\begin{equation}
  \begin{split}
    Q(\boldsymbol\beta, \boldsymbol\beta^g) & = \sum_{l = 1}^M \sum_{i = 1}^N \ln (\alpha_l p_l (\mathbf{x}_i \mid \theta_l)) p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) \\
    & = \sum_{l = 1}^M \sum_{i = 1}^N (\ln \alpha_l) p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) + \sum_{l = 1}^M \sum_{i = 1}^N (\ln p_l (\mathbf{x}_i \mid \theta_l)) p(l \mid \mathbf{x}_i, \boldsymbol\beta^g)
  \end{split}
  (\#eq:mixestep)
\end{equation}

Now M-step: maximize $Q$.

\begin{equation}
  \begin{cases}
    \alpha_l^{new} = \frac{1}{N} \sum_{i = 1}^N p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) \\
    \boldsymbol\mu_l^{new} = \frac{\sum \mathbf{x}_i p(l \mid \mathbf{x}_i, \boldsymbol\beta^g)}{\sum p(l \mid \mathbf{x}_i, \boldsymbol\beta^g)} \in \R^p \\
    \Sigma_l^{new} = \frac{\sum p(l \mid \mathbf{x}_i, \boldsymbol\beta^g) (\mathbf{x}_l - \boldsymbol\mu_l^{new}) (\mathbf{x}_l - \boldsymbol\mu_l^{new})^T}{\sum p(l \mid \mathbf{x}_i, \boldsymbol\beta^g)} \in \R^{p \times p}
  \end{cases}
  (\#eq:mixmstep)
\end{equation}

Refer to Example \@ref(exm:mixgam). Use the same generated data.

```{r}
mix_em <- function(fn = dgamma, x, N = 10000, par, par_name = "rate", tol = .Machine$double.eps^.5, ...) {
  ll <- list(x = x, ...)
  dens <- 
    tibble(
      key = paste0("f", 1:3),
      mu = double(3L),
      lam = par
    )
  lam_old <- par + 1
  for (i in seq_len(N)) {
    dens <- 
      lapply(
        dens$lam,
        function(y) {
          ll[[par_name]] <- 1 / (2 * y)
          do.call(fn, ll)
        }
      ) %>% 
      bind_cols() %>% 
      rename_all(
        .funs = list(
          ~str_replace_all(., pattern = "V", replacement = "f")
        )
      ) %>% 
      mutate(total = apply(., 1, sum)) %>% 
      gather(-total, key = "key", value = "value") %>% 
      mutate(post = value / total) %>% # posterior prob y from each f1, f2, f3 - E step
      group_by(key) %>% 
      summarise(mu = sum(y * post) / sum(post)) %>%  # update means - M step
      mutate(lam = mu / sum(mu))
    if (
      dens %>%
      summarise(
        sum(abs(lam - lam_old) / lam_old)
      ) %>%
      pull() < tol
    ) break()
    lam_old <- dens %>% select(lam) %>% pull()
  }
  list(lambda = dens %>% select(lam), iter = i)
}
```

The result for same initial $\lambda$ is:

```{r}
mix_em(x = y, par = c(.5, .4, .1), shape = 1 / 2)
```