+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+01_ Mixture 3 Gaussian_1 Dimention
+Tran Le
+2/24/2022
+ +
+
+1 Introduction
+Here is the R code that I use to illustrate a simulation example using the EM algorithm for Gaussian mixture models with 3 (Gaussian) distributions in 1D. The EM formula construction and the pseudo-code can be found here.
+
+
+2 Utility functions
+-
+
- Functions to create and plot the synthesis data. +
creat_3G_dat_f <- function(n, p1, mu1, sig1, p2, mu2, sig2, mu3, sig3){
+ # # Create synthesis data
+ # n: the total sample size of the three distributions
+ # p1, p2: the probability of a data point belongs to the 1st, 2nd distribution
+ # p1+p2 < 1 since we have p3=1-(p1+p2) is the proportion of 3rd distribution
+ n1 <- floor(n*p1); Y1 <- rnorm(n1, mu1, sig1)
+ n2 <- floor(n*p2); Y2 <- rnorm(n2, mu2, sig2)
+ n3 <- ceiling(n*(1-p1-p2)); Y3 <- rnorm(n3, mu3, sig3)
+ data <- data.frame(
+ dist = c( rep("1", n1), rep("2", n2) , rep("3", n3)),
+ value = c( Y1, Y2, Y3))
+ return(data)
+}
+
+plot_3G_syn_dat_f <- function(data){
+ # a data set created by function creat_dat_f
+ p1 <- data %>%
+ ggplot(aes(x=value, fill=dist), addDensityCurve=TRUE) +
+ geom_histogram( color="#e9ecef", alpha=0.6,
+ position = 'identity', binwidth=1) +
+ scale_fill_manual(values=c("#69b3a2", "#404080", "#80427d"))
+ p2 <- data %>%
+ ggplot(aes(x=value)) +
+ geom_histogram( color="#e9ecef", alpha=0.6, position = 'identity', binwidth=1)
+ theplot <- plot_grid(p1, p2, labels=c("Three data components","Mixture data"),
+ ncol = 2, nrow = 1)
+ return(theplot)
+}-
+
- Function to update the posterior probability +
com_3G_pi_f <-function(y, p1, mu1, sig1, p2, mu2, sig2, mu3, sig3){
+ # compute posterior probability associative with each yi
+ f1 <- dnorm(y, mean = mu1, sd = sig1)
+ f2 <- dnorm(y, mean = mu2, sd = sig2)
+ f3 <- dnorm(y, mean = mu3, sd = sig3)
+ denorm <- (p1*f1)+(p2*f2)+((1-p1-p2)*f3)
+ pi_i1 <- (p1*f1)/denorm
+ pi_i2 <- (p2*f2)/denorm
+ return(c(pi_i1, pi_i2, 1-pi_i1-pi_i2))
+}-
+
- The main EM function +
EM_Three_gaussian_f <-function(Y, p1, mu1, sig1, p2, mu2, sig2, mu3, sig3, eps, max_iter, verbose) {
+ # Function to estimate parameters based on initial guesses
+ # p1, mu1, sig1, p2, mu2, sig2, mu3, sig3: INITIAL parameter values
+ # eps" stopping criterion threshold
+ # max_iter: maximum number of EM iteration
+ # verbose = T if you want to see the updated estimation values in each iteration step
+ n <- length(Y)
+ for (i in (1:max_iter)){
+ # Compute the (s+1) step from (s) step
+ pij <- map(Y, ~com_3G_pi_f(.x, p1, mu1, sig1, p2, mu2, sig2, mu3, sig3))
+ pij <- as.data.frame(do.call(rbind, pij))
+ p1_new <- sum(pij$V1)/n
+ p2_new <- sum(pij$V2)/n
+ #
+ mu1_new <- sum(Y*pij$V1)/sum(pij$V1)
+ sig1_new <- sqrt(sum(((Y-mu1)^2)*pij$V1)/sum(pij$V1))
+ mu2_new <- sum(Y*pij$V2)/sum(pij$V2)
+ sig2_new <- sqrt(sum(((Y-mu2)^2)*pij$V2)/sum(pij$V2))
+ mu3_new <- sum(Y*pij$V3)/sum(pij$V3)
+ sig3_new <- sqrt(sum(((Y-mu3)^2)*pij$V3)/sum(pij$V3))
+ if(verbose){cat("\n at step ", i, ": ",
+ mu1_new, sig1_new, mu2_new,sig2_new, mu3_new, sig3_new)}
+ # stopping criteria
+ if (max(abs(p1_new-p1), abs(p2_new -p2),
+ abs(mu1_new -mu1), abs(mu2_new-mu2), abs(mu3_new-mu3))< eps){break}
+ # Update
+ p1 <- p1_new; mu1 <- mu1_new; sig1 <- sig1_new
+ p2 <- p2_new; mu2 <- mu2_new; sig2 <- sig2_new
+ mu3 <- mu3_new; sig3 <- sig3_new
+ } # end for loop
+ if (i==max_iter){cat("\n The stopping criterium does not meet \n")
+ } else {cat("\n The EM converge at iteration", i, '\n')}
+ return(list(pij=pij,
+ p1=p1_new, mu1=mu1_new, sig1=sig1_new,
+ p2=p2_new, mu2=mu2_new, sig2=sig2_new,
+ p3= 1- p1_new-p2_new, mu3=mu3_new, sig3=sig3_new))
+}-
+
- A function to do the EM clustering. +
create_comp_dat_3G_1D_f <- function(data, pij){
+ # Out: comp_dat_2D: contain original label (dist) and EM cluster prediction (pred)
+ pij$Max<-colnames(pij)[apply(pij,1,which.max)]
+ pij <- pij %>% mutate(Max_col = str_sub(Max, start = 2, end = 2))
+ comp_dat_2D <- cbind(data, pred = pij$Max_col)
+ comp_dat_2D <- comp_dat_2D %>% mutate(err = as.integer(dist)-as.integer(pred))
+ return(comp_dat_2D)
+}-
+
- Function to plot the original data and the predicted data (histograms) by EM clustering. +
plot_1D_org_pre_dat_f <- function(comp_dat){
+ # Plot the original and predicted data gotten from EM algorithm
+ p1 <- comp_dat %>%
+ ggplot(aes(x=value, fill=dist), addDensityCurve=TRUE) +
+ geom_histogram( color="#e9ecef", alpha=0.6,
+ position = 'identity', binwidth=1) +
+ scale_fill_manual(values=c("#69b3a2", "#404080", "#80427d"))
+ p2 <- comp_dat %>%
+ ggplot(aes(x=value, fill=pred)) +
+ geom_histogram( color="#e9ecef", alpha=0.6,
+ position = 'identity', binwidth=1)
+ theplot <- plot_grid(p1, p2,
+ labels=c("Orginial data components",
+ "Predicted data components"),
+ ncol = 2, nrow = 1)
+ return(theplot)
+}
+
+
+
+
+3 Example
+# EXAMPLE 1 Set up a model----------------------------------------
+set.seed(538725)
+# Parameters true values
+Tv <- data.frame(p1=0.5, mu1=50, sig1=3,
+ p2=0.2, mu2=80, sig2=4,
+ mu3=85, sig3=3)
+data <- creat_3G_dat_f(n=1000,
+ p1=Tv$p1, mu1=Tv$mu1, sig1=Tv$sig1,
+ p2=Tv$p2, mu2=Tv$mu2, sig2=Tv$sig2,
+ mu3=Tv$mu3, sig3=Tv$sig3)
+# png("./EX_Result/EX1_3G_1d_data.png");
+plot_3G_syn_dat_f(data)
# ;dev.off()
+# Parameter Initial values
+
+3.1 If we choose the correct number of clusters
+Iv <- data.frame(p1=0.3, mu1=50, sig1=sd(data$value)/3,
+ p2=0.3, mu2=75, sig2=sd(data$value)/3,
+ mu3=89, sig3=sd(data$value)/3)
+# Solve the problem using Em
+EX3 <- EM_Three_gaussian_f(Y=data$value,
+ p1=Iv$p1, mu1=Iv$mu1, sig1=sd(data$value)/3,
+ p2=Iv$p2, mu2=Iv$mu2, sig2=sd(data$value)/3,
+ mu3=Iv$mu3, sig3=sd(data$value)/3,
+ eps=0.02, max_iter=1000, verbose=F) ##
+## The EM converge at iteration 20# Get the comparison data
+EX3_pij <- as.data.frame(EX3$pij)
+EX3_comp_dat <- create_comp_dat_3G_1D_f(data, EX3_pij)
+# png("./EX_Result/EX1_3G_1d_orgpre.png");
+plot_1D_org_pre_dat_f(EX3_comp_dat)
# ;dev.off()
+EX3_error <- tabyl(EX3_comp_dat$err) %>% mutate(error = `EX3_comp_dat$err`) %>%
+ select(c("error", "n", "percent"))
+# get the parameter estimator
+EX3 <- as.data.frame(EX3[c(2:10)])
+save(Tv, Iv, EX3, EX3_error, file= "./EX_Result/EX3.RData")
+
+3.2 If we choose a wrong number of clusters
+# If we use the #clusters =2
+source("01_Mixture_2Gaussian_1d.R")
+EX3_2 <- EM_two_gaussian_f(mu1_0 = 50, sig1_0= 10,
+ mu2_0=85, sig2_0=10,
+ p_0=0.5, Y= data$value,
+ max_iter=1000,
+ eps=0.01, dist_thres_prob= 0.5, verbose=F) %>%
+ as.data.frame()
+U_vec <- com_U_vec_f(dist_thres_prob=0.5, EX3_2$p, Y=data$value,
+ mu1=EX3_2$mu1, sig1=EX3_2$sig1,
+ mu2=EX3_2$mu2, sig2=EX3_2$sig2)
+EX3_2_comdat <- cbind(data, U_vec) %>%
+ mutate(dist = as.numeric(dist)) %>%
+ mutate(err= dist-U_vec)
+EX3_2_comdat_err <- tabyl(EX3_2_comdat$err)
+EX3_2## i p mu1 sig1 mu2 sig2
+## 1 4 0.5 49.94347 3.073976 82.99561 4.113227EX3_2_comdat_err## EX3_2_comdat$err n percent
+## 0 700 0.7
+## 1 300 0.3# write_rds(EX3_2, "./EX_Result/EX3_2.rds")
+# write_rds(EX3_2_comdat_err, "./EX_Result/EX3_2_comdat_err.rds")plot_org_pre_wrongk_dat_f <- function(comp_dat){
+ # Plot the original + predicted data gotten from EM algorithm
+ # when we choose a WRONG number of cluster
+ comp_dat <- comp_dat %>% mutate(across(c("dist", "U_vec"), as.factor)) %>%
+ rename(pre_dist = U_vec)
+ p1 <- comp_dat %>%
+ ggplot(aes(x=value, fill=dist), addDensityCurve=TRUE) +
+ geom_histogram( color="#e9ecef", alpha=0.6, position = 'identity', binwidth=1) +
+ scale_fill_manual(values=c("#69b3a2", "#404080", "#80427d"))
+ p2 <- comp_dat %>%
+ ggplot(aes(x=value, fill=pre_dist), addDensityCurve=TRUE) +
+ geom_histogram( color="#e9ecef", alpha=0.6, position = 'identity', binwidth=1) +
+ scale_fill_manual(values=c("#69b3a2", "#404080"))
+ theplot <- plot_grid(p1, p2,
+ labels=c("Orginial data",
+ "Predicted data"), ncol = 2, nrow = 1)
+ return(theplot)
+}
+
+# png("./EX_Result/EX3_wrong_k.png");
+plot_org_pre_wrongk_dat_f(EX3_2_comdat)
# ;dev.off()
+
+