Gaussian mixture

Author

Affiliation

Beniamino Sartini

University of Bologna

library(dplyr)
# required for figures  
library(ggplot2)
library(gridExtra)
# required for render latex 
library(backports)
library(latex2exp)
# required for render tables
library(knitr)
library(kableExtra)
# Random seed 
set.seed(1)

Let’s consider a linear combination of a Bernoulli and two normal random variables, all assumed to be independent, i.e. $$\begin{array}{}\text{(1)}& X\sim B\cdot Z_1+(1-B)\cdot Z_2\text{,}\end{array}$$ or in compact form $X\sim GM({\mu }_1,{\mu }_2,{\sigma }_1^2,{\sigma }_2^2,p)$. Formally, $B\sim \text{Ber}(p)$ is a Bernoulli while $Z_1=\mathcal{N}({\mu }_1,{\sigma }_1^2)$ and $Z_2=\mathcal{N}({\mu }_2,{\sigma }_2^2)$ are two independent Gaussian random variables.

# ================== Setups ==================
t_bar <- 5000 # number of steps ahead  
# parameters
par <- c(mu1=-2,mu2=2,sd1=1,sd2=1,p=0.5) 
# ============================================
# Gaussian mixture simulation 
N_1 <- rnorm(t_bar, mean = par[1], sd = par[3])
N_2 <- rnorm(t_bar, mean = par[2], sd = par[4])
B <- rbinom(t_bar, 1, prob = par[5])
Xt <- B*N_1 + (1 - B)*N_2
# Empiric pdf and cdf  
ker <- density(Xt, from = min(Xt), to = max(Xt))
ker$cdf_emp <- cumsum(ker$y/sum(ker$y))
# Components normal pdf 
ker$pdf_Z1 <- dnorm(ker$x, mean = par[1], sd = par[3]) 
ker$pdf_Z2 <- dnorm(ker$x, mean = par[2], sd = par[4])
# Mixture pdf and cdf 
ker$pdf <- par[5]*ker$pdf_Z1 + (1-par[5])*ker$pdf_Z2
ker$cdf <- cumsum(ker$pdf/sum(ker$pdf))
# =================== Plot ===================
# Plot trajectory
plot_gm <- ggplot()+
  geom_point(aes(1:t_bar, Xt), alpha = exp(-0.00009*t_bar))+
  labs(x = "t", y = TeX("$X_t$"))+
  theme_bw()
# Plot pdf 
plot_pdf <- ggplot()+
  geom_line(aes(ker$x, ker$y))+
  geom_line(aes(ker$x, ker$pdf), color = "red")+
  labs(x = NULL, y = "Pdf")+
  theme_bw()+
  coord_flip()
# Plot cdf 
plot_cdf <- ggplot()+
  geom_line(aes(ker$x, ker$cdf_emp))+
  geom_line(aes(ker$x, ker$cdf), color = "red")+
  labs(x = NULL, y = "Cdf")+
  theme_bw()+
  coord_flip()
plot_gm
gridExtra::grid.arrange(plot_pdf, plot_cdf, ncol = 2)

(a) Simulation.

(b) Pdf and Cdf.

Figure 1: Gaussian Mixture simulation and density function with parameters ${\mu }_1=-2$, ${\mu }_2=2$, ${\sigma }_1=1$, ${\sigma }_2=1$ and $p=0.5$.

1 Distribution and density

The distribution function of a Gaussian mixture reads explicitely as: $$\begin{array}{}\text{(2)}& F_X(x)=p\cdot \mathrm{\Phi }\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+(1-p)\cdot \mathrm{\Phi }\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\text{,}\end{array}$$ where $\mathrm{\Phi }$ is the cumulative distribution function of a standard normal random variable. Taking the derivative, it can be easily shown that the density function reads: $$\begin{array}{}\text{(3)}& f_X(x)=p\cdot \varphi \left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+(1-p)\cdot \varphi \left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\text{.}\end{array}$$ where $\varphi$ is the density function of a standard normal random variable.

R

dnorm_mix <- function(params) {
  # parameters
  mu1 = params[1]
  mu2 = params[2]
  sd1 = params[3]
  sd2 = params[4]
  p = params[5]
  function(x, log = FALSE){
    probs <- p*stats::dnorm(x, mean = mu1, sd = sd1) + (1-p)*stats::dnorm(x, mean = mu2, sd = sd2)
    if (log) {
      probs <- base::log(probs)
    }
    return(probs)
  }
}

Density of a Gaussian Mixture

Proof. The distribution function of a Gaussian Mixture is defined as: $$F_X(y)=\mathbb{P}(X\le y)=\mathbb{E}\{{\mathbb{1}}_{X\le x}\}$$ Hence, we can rewrite it in terms of conditional expectation with respect to $B$, i.e.  $$\begin{array}{rl}{F}_X(y)& =\mathbb{E}\{{\mathbb{1}}_{X\le x}|B\}=\\ & =\mathbb{E}\{{\mathbb{1}}_{X\le x}|B=0\}\mathbb{P}(B=0)+\mathbb{E}\{{\mathbb{1}}_{X\le x}|B=1\}\mathbb{P}(B=1)=\\ & =p\cdot \mathbb{P}(Z_1\le x)+(1-p)\cdot \mathbb{P}(Z_2\le x)\end{array}$$ Hence, standardizing the normal random variable we obtain, $$F_X(x)=p\cdot \mathrm{\Phi }\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+(1-p)\cdot \mathrm{\Phi }\left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\text{.}$$ where $\mathrm{\Phi }$ denotes the distribution function of a standard normal. Knowing that $f_X(x)=\frac{d{F}_X(x)}{dx}$ and that ${\varphi }_X(x)=\frac{d\mathrm{\Phi }(x)}{dx}$, where $\varphi$ is the density function of a standard normal we obtain the result, i.e. $$f_X(x)=p\cdot \varphi \left(\frac{x-{\mu }_1}{{\sigma }_1}\right)+(1-p)\cdot \varphi \left(\frac{x-{\mu }_2}{{\sigma }_2}\right)\text{.}$$

2 Moments

Given that $Z_1$, $Z_2$ and $B$ are independent, the expectation can be computed as: $$\mathbb{E}\{X\}=p{\mu }_1+(1-p){\mu }_2$$ The second moment is computed as: $$\mathbb{E}\{X^2\}=p({\mu }_1^2+{\sigma }_1^2)+(1-p)({\mu }_2^2+{\sigma }_2^2)$$ Hence, the variance reads: $$\begin{array}{r}\mathbb{V}\{X\}=p(1-p)({\mu }_1-{\mu }_2{)}^2+{\sigma }_1^{2}p+{\sigma }_2^2(1-p)\end{array}$$

Gaussian mixture moments

Proof. Given that $Z_1$, $Z_2$ and $B$ are independent, the expectation can be computed as: $$\begin{array}{rl}\mathbb{E}\{X\}& =\mathbb{E}\{\mathbb{E}\{X|B\}\}=\\ & =\mathbb{E}\{X|B=1\}\mathbb{P}(B=1)+\mathbb{E}\{X|B=0\}\mathbb{P}(B=0)=\\ & =p\mathbb{E}\{Z_1\}+(1-p)\mathbb{E}\{Z_2\}=\\ & =p{\mu }_1+(1-p){\mu }_2\end{array}$$ The second moment is computed similarly to the first one, i.e.  $$\begin{array}{rl}\mathbb{E}\{X^2\}& =\mathbb{E}\{\mathbb{E}\{X^2|B\}\}=\\ & =\mathbb{E}\{X^2|B=1\}\mathbb{P}(B=1)+\mathbb{E}\{X^2|B=0\}\mathbb{P}(B=0)=\\ & =\mathbb{E}\{B\}\mathbb{E}\{Z_1^2\}+\mathbb{E}\{1-B\}\mathbb{E}\{Z_2^2\}=\\ & =p\mathbb{E}\{Z_1^2\}+(1-p)\mathbb{E}\{Z_2^2\}=\\ & =p({\mu }_1^2+{\sigma }_1^2)+(1-p)({\mu }_2^2+{\sigma }_2^2)\end{array}$$ The variance, by definition, is given by: $$\mathbb{V}\{X\}=\mathbb{E}\{X^2\}-\mathbb{E}\{X{\}}^2$$ where the first moment squared is $$\begin{array}{rl}\mathbb{E}\{X{\}}^2& ={[p{\mu }_1+(1-p){\mu }_2]}^2=\\ & =p^2{\mu }_1^2+(1-p{)}^2{\mu }_2^2+2p(1-p){\mu }_1{\mu }_2\end{array}$$ Hence the variance, $$\begin{array}{rl}\mathbb{V}\{X\}& =p({\mu }_1^2+{\sigma }_1^2)+(1-p)({\mu }_2^2+{\sigma }_2^2)-p^2{\mu }_1^2-(1-p{)}^2{\mu }_2^2-2p(1-p){\mu }_1{\mu }_2=\\ & =p{\mu }_1^2+p{\sigma }_1^2+{\mu }_2^2+{\sigma }_2^2-p{\mu }_2^2-p{\sigma }_2^2-p^2{\mu }_1^2-(1-p{)}^2{\mu }_2^2-2p(1-p){\mu }_1{\mu }_2=\\ & ={\mu }_1^{2}p(1-p)+p{\sigma }_1^2+(1-p){\sigma }_2^2+p(1-p){\mu }_2^2-2p(1-p){\mu }_1{\mu }_2=\\ & =p(1-p)({\mu }_1^2-{\mu }_2^2-2{\mu }_1{\mu }_2)+p{\sigma }_1^2+(1-p){\sigma }_2^2=\\ & =p(1-p)({\mu }_1-{\mu }_2{)}^2+p{\sigma }_1^2+(1-p){\sigma }_2^2\end{array}$$

3 Maximum likelihood

Minimizing the negative log-likelihood gives an estimate of the parameters, i.e. $$\underset{{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,p}{\text{argmin}}\{\sum _{i=1}^t\log (f_X(x_i))\}\text{,}$$ or equivalently maximizing the negative log-likelihood, i.e.  $$\underset{{\mu }_1,{\mu }_2,{\sigma }_1,{\sigma }_2,p}{\text{argmax}}\{-\sum _{i=1}^t\log (f_X(x_i))\}\text{.}$$

Example: ML-estimate

# Initialize parameters 
init_params <- par*runif(5, 0.3, 1.1)
# Log-likelihood function
log_lik <- function(params, x){
  # Parameters
  mu1 = params[1]
  mu2 = params[2]
  sd1 = params[3]
  sd2 = params[4]
  p = params[5]
  # Ensure that probability is in (0,1)
  if(p > 0.99 | p < 0.01 | sd1 < 0 | sd2 < 0){ s
    return(NA_integer_)
  }
  # Mixture density
  pdf_mix <- dnorm_mix(params)
  # Log-likelihood
  loss <- sum(pdf_mix(x, log = TRUE), na.rm = TRUE)
  return(loss)
}
# Optimal parameters
# fnscale = -1 to maximize (or use negative likelihood)
ml_estimate <- optim(par = init_params, log_lik, 
                     x = Xt, control = list(maxit = 500000, fnscale = -1))

Table 1: Maximum likelihood estimates for a Gaussian Mixture.

Parameter	True	Estimate	Log-lik	Bias
$${\mu }_1$$	-2.0	-1.9905803	-10332.56	0.0094197
$${\mu }_2$$	2.0	2.0006381	-10332.56	0.0006381
$${\sigma }_1$$	1.0	1.0457653	-10332.56	0.0457653
$${\sigma }_2$$	1.0	1.0033373	-10332.56	0.0033373
$$p$$	0.5	0.4937459	-10332.56	-0.0062541

4 Moments matching

Let’s fix the parameter of the first component, namely ${\mu }_1$ and ${\sigma }_1^2$ and a certain probability $p$. Then let’s compute the sample estimate of the expectation of $X_t$, i.e. $$\mathbb{E}\{X\}=\frac{1}{t}\sum _{i=1}^t{x}_i=\hat{\mu }$$ and the sample variance: $$\mathbb{V}\{X\}=\frac{1}{t}\sum _{i=1}^t(x_i-\hat{\mu }{)}^2={\hat{\sigma }}^2$$ In order to obtain an estimate of the second distribution such that the Gaussian Mixture moments match exactly the sample estimates we solve the system for ${\mu }_2$ and ${\sigma }_2^2$: $$\{\begin{array}{l}\hat{\mu }=p{\mu }_1+(1-p){\mu }_2\\ {\hat{\sigma }}^2=p(1-p)({\mu }_1-{\mu }_2{)}^2+{\sigma }_1^{2}p+{\sigma }_2^2(1-p)\end{array}$$ which lead to a unique solution, i.e. $$\begin{array}{rl}& {\mu }_2=\frac{\hat{\mu }-p{\mu }_1}{1-p}\\ & {\sigma }_2^2=\frac{{\hat{\sigma }}^2-p{\sigma }_1^2}{1-p}-p({\mu }_1-{\mu }_2{)}^2\end{array}$$

5 Moment generating function

The moment generating function of a Gaussian mixture random variable (Equation 1) reads: $$M_X(t)=p\cdot \exp \{{\mu }_{1}t+\frac{t^2{\sigma }_1^2}{2}\}+(1-p)\cdot \exp \{{\mu }_{2}t+\frac{t^2{\sigma }_2^2}{2}\}$$

Moment generating function

Proof. Applying the definition of moment generating function on a gaussian mixture we obtain: $$\begin{array}{rl}{M}_X(t)& =\mathbb{E}\{e^{tX}\}=\\ & =p\cdot \mathbb{E}\{e^{t{Z}_1}\}+(1-p)\cdot \mathbb{E}\{e^{t{Z}_2}\}=\\ & =p\cdot M_{Z_1}(t)+(1-p)\cdot M_{Z_1}(t)=\\ & =p\cdot \exp \{{\mu }_{1}t+\frac{t^2{\sigma }_1^2}{2}\}+(1-p)\cdot \exp \{{\mu }_{2}t+\frac{t^2{\sigma }_2^2}{2}\}\end{array}$$ where $M_{Z_1}(t)$ and $M_{Z_2}(t)$ are the moment generating functions of a Gaussian random variable.

6 Esscher transform

The Esscher transform of a Gaussian mixture random variable reads: $$E_{\theta }\{f_X(x)\}=\tilde{p}\cdot \mathcal{N}({\mu }_1+\theta {\sigma }_1^2,{\sigma }_1^2)+(1-\tilde{p})\cdot \mathcal{N}({\mu }_2+\theta {\sigma }_2^2,{\sigma }_2^2)$$ where $\tilde{p}$ are defined as: $$\tilde{p}=\frac{p{M}_{Z_1}(\theta )}{p{M}_{Z_1}(\theta )+(1-p)M_{Z_2}(\theta )}$$

Moment generating function

Proof. In general, the Esscher transform of a density function is computed as:
$$E_{\theta }\{f_X(x)\}=\frac{e^{\theta X}{f}_X(x)}{M_X(\theta )}=\frac{e^{\theta X}{f}_X(x)}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\theta X}{f}_X(x)dx}$$ Substituing the density function of a Gaussian mixture we obtain: $$E_{\theta }\{f_X(x)\}=\frac{e^{\theta X}(p{f}_{Z_1}(x)+(1-p)f_{Z_2}(x))}{{\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{e}^{\theta X}(p{f}_{Z_1}(x)+(1-p)f_{Z_2}(x))dx}$$ Let’s examine the product of the first component: $$p{e}^{\theta X}{f}_{Z_1}(x)=\frac{p{e}^{\theta x}}{\sqrt{2\pi }{\sigma }_1}\exp \{-\frac{(x-{\mu }_1{)}^2}{2{\sigma }_1^2}\}=\frac{1}{\sqrt{2\pi }{\sigma }_1^2}\exp \{-\frac{(x-{\mu }_1{)}^2}{2{\sigma }_1^2}+\theta x\}$$ In order to rewrite it as the pdf of a normal with mean $({\mu }_1+\theta {\sigma }_1^2)$, we note that the expansion of such product reads: $$(x-{\mu }_1-\theta {\sigma }_1^2{)}^2=x^2-2x{\mu }_1-2\theta x{\sigma }_1^2+{\mu }^2+2\mu \theta {\sigma }_1^2+{\theta }^2{\sigma }_1^2$$ Hence, let’s add and subtract inside the exponential $\pm {\mu }_1\theta \pm \frac{{\theta }^2{\sigma }_1^2}{2}$, i.e. $$\begin{array}{rl}p{e}^{\theta X}{f}_{Z_1}(x)& =\frac{p}{\sqrt{2\pi }{\sigma }_1}\exp \{-\frac{(x-{\mu }_1{)}^2}{2{\sigma }_1^2}+\theta x-{\mu }_1\theta -\frac{{\theta }^2{\sigma }_1^2}{2}\}\exp \{{\mu }_1\theta +\frac{{\theta }^2{\sigma }_1^2}{2}\}=\\ & =\frac{p}{\sqrt{2\pi }{\sigma }_1}\exp \{{\mu }_1\theta +\frac{\theta {\sigma }_1^2}{2}\}\exp \{-\frac{(x-{\mu }_1-{\theta }^2{\sigma }_1^2{)}^2}{2{\sigma }_1^2}\}\end{array}$$ Hence let’s collect the first part and the denominator (mgf of the Gaussian mixture in $\theta$) and define the new probability: $$\tilde{p}=\frac{p\exp \{{\mu }_1\theta +\frac{{\theta }^2{\sigma }_1^2}{2}\}}{p\exp \{{\mu }_1\theta +\frac{{\theta }^2{\sigma }_1^2}{2}\}+(1-p)\exp \{{\mu }_2\theta +\frac{{\theta }^2{\sigma }_2^2}{2}\}}$$ That can be equivalently written in terms of the moment generating functions, i.e.  $$\tilde{p}=\frac{p{M}_{Z_1}(\theta )}{p{M}_{Z_1}(\theta )+(1-p)M_{Z_2}(\theta )}$$ Hence, repeating the same for the second component we obtain that the Esscher transform of a Gaussian mixture is defined as: $$E_{\theta }\{f_X(x)\}=\tilde{p}\cdot \mathcal{N}({\mu }_1+\theta {\sigma }_1^2,{\sigma }_1^2)+(1-\tilde{p})\cdot \mathcal{N}({\mu }_2+\theta {\sigma }_2^2,{\sigma }_2^2)$$