Question

For expert using R

I try to solve this question((USING DATA FAITHFUL)) but each time I solve it, I have error , I try it many times. So,everything you write will be helpful..

Modify the EM-algorithm functions to work for a general K component Gaussian mixtures. Please use this function to fit a K= 1;2;3;4 modelto the old faithful data available in R (You need to initialize the EM-algorithm First ).   

Which modelseems to t the data better? (Hint: use BIC to compare models.)

Here what I try to use

## EM algorithm for univariate normal mixture

# The E-step

E.step <- function(x, pi, Mu, S2){

K <- length(pi)

n <- length(x)

tau <- matrix(rep(NA, n * K), ncol = K)

for (i in 1:n){

for (k in 1:K){

tau[i,k] <- pi[k] * dnorm(x[i], Mu[k], sqrt(S2[k]))

}

tau[i,] <- tau[i,] / sum(tau[i,])

}

return(tau)

}

#The M-step

M.step <- function(x, tau){

n <- length(x)

K <- dim(tau)[2]

tau.sum <- apply(tau, 2, sum)

pi <- tau.sum / n

Mu <- t(tau) %*% x / tau.sum

S2[1] <- t(tau[,1]) %*% (x - Mu[1])^2 / tau.sum[1]

S2[2] <- t(tau[,2]) %*% (x - Mu[2])^2 / tau.sum[2]

return(list(pi = pi, Mu = Mu, S2 = S2))

}

## The log-likelihood function

logL <- function(x, pi, Mu, S2){

n <- length(x)

ll <- 0

for (i in 1:n){

ll <- ll + log(pi[1] * dnorm(x[i], Mu[1], sqrt(S2[1])) +

pi[2] * dnorm(x[i], Mu[2], sqrt(S2[2])))

}

return(ll)

}

## The algorithm

EM <- function(x, pi, Mu, S2, tol){

t <- 0

ll.old <- -Inf

ll <- logL(x, pi, Mu, S2)

repeat{

t <- t + 1

if ((ll - ll.old) / abs(ll) < tol) break

ll.old <- ll

tau <- E.step(x, pi, Mu, S2)

M <- M.step(x, tau)

pi <- M$pi

Mu <- M$Mu

S2 <- M$S2

ll <- logL(x, M$pi, M$Mu, M$S2)

cat("Iteration", t, "logL =", ll, " ")

}

return(list(pi = M$pi, Mu = M$Mu, S2 = M$S2, tau = tau, logL = ll))

}

## generate data

set.seed(1)

pi <- c(0.3, 0.7)

Mu <- c(5, 10)

S2 <- c(1, 1)

n <- 1000

n1 <- rbinom(1, n, pi[1])

n2 <- n - n1

x1 <- rnorm(n1, Mu[1], sqrt(S2[1]))

x2 <- rnorm(n2, Mu[2], sqrt(S2[2]))

x <- c(x1, x2)

hist(x, freq = FALSE, ylim = c(0, 0.2))

# pick initial values

pi.init <- c(0.5, 0.5)

Mu.init <- c(3, 10)

S2.init <- c(0.4, 2)

#Run EM

A <- EM(x, pi.init, Mu.init, S2.init, tol = 10^-6)

#plot

t <- seq(0, 15, by = 0.01)

y <- pi[1] * dnorm(t, Mu[1], sqrt(S2[1])) +

pi[2] * dnorm(t, Mu[2], sqrt(S2[2]))

y.est <- A$pi[1] * dnorm(t, A$Mu[1], sqrt(A$S2[1])) +

A$pi[2] * dnorm(t, A$Mu[2], sqrt(A$S2[2]))

points(t, y, type = "l")

points(t, y.est, type = "l", col = 2, lty = 2)

# assign observations to components - clustering

d <- function(x) which(x == max(x))

apply(A$tau, 1, d)

apply(A$tau, 1, which.max)

# assess misclassification

table(apply(A$tau, 1, which.max), c(rep(1, n1), rep(2, n2)))

Answer 1

there is nothing wrong in your code. you just missed a space in the last line of the code.

table(apply(A $tau, 1, which.max), c(rep(1, n1), rep(2, n2)))

make sure that you use R programming version 3.5.1

the table obtained at the last step is

	1	2
1	296	3
2	5	696

execute with this code

## EM algorithm for univariate normal mixture

# The E-step

E.step <- function(x, pi, Mu, S2){

K <- length(pi)

n <- length(x)

tau <- matrix(rep(NA, n * K), ncol = K)

for (i in 1:n){

for (k in 1:K){

tau[i,k] <- pi[k] * dnorm(x[i], Mu[k], sqrt(S2[k]))

}

tau[i,] <- tau[i,] / sum(tau[i,])

}

return(tau)

}

#The M-step

M.step <- function(x, tau){

n <- length(x)

K <- dim(tau)[2]

tau.sum <- apply(tau, 2, sum)

pi <- tau.sum / n

Mu <- t(tau) %*% x / tau.sum

S2[1] <- t(tau[,1]) %*% (x - Mu[1])^2 / tau.sum[1]

S2[2] <- t(tau[,2]) %*% (x - Mu[2])^2 / tau.sum[2]

return(list(pi = pi, Mu = Mu, S2 = S2))

}

## The log-likelihood function

logL <- function(x, pi, Mu, S2){

n <- length(x)

ll <- 0

for (i in 1:n){

ll <- ll + log(pi[1] * dnorm(x[i], Mu[1], sqrt(S2[1])) +

pi[2] * dnorm(x[i], Mu[2], sqrt(S2[2])))

}

return(ll)

}

## The algorithm

EM <- function(x, pi, Mu, S2, tol){

t <- 0

ll.old <- -Inf

ll <- logL(x, pi, Mu, S2)

repeat{

t <- t + 1

if ((ll - ll.old) / abs(ll) < tol) break

ll.old <- ll

tau <- E.step(x, pi, Mu, S2)

M <- M.step(x, tau)

pi <- M$pi

Mu <- M$Mu

S2 <- M$S2

ll <- logL(x, M$pi, M$Mu, M$S2)

cat("Iteration", t, "logL =", ll, " ")

}

return(list(pi = M$pi, Mu = M$Mu, S2 = M$S2, tau = tau, logL = ll))

}

## generate data

set.seed(1)

pi <- c(0.3, 0.7)

Mu <- c(5, 10)

S2 <- c(1, 1)

n <- 1000

n1 <- rbinom(1, n, pi[1])

n2 <- n - n1

x1 <- rnorm(n1, Mu[1], sqrt(S2[1]))

x2 <- rnorm(n2, Mu[2], sqrt(S2[2]))

x <- c(x1, x2)

hist(x, freq = FALSE, ylim = c(0, 0.2))

# pick initial values

pi.init <- c(0.5, 0.5)

Mu.init <- c(3, 10)

S2.init <- c(0.4, 2)

#Run EM

A <- EM(x, pi.init, Mu.init, S2.init, tol = 10^-6)

#plot

t <- seq(0, 15, by = 0.01)

y <- pi[1] * dnorm(t, Mu[1], sqrt(S2[1])) +

pi[2] * dnorm(t, Mu[2], sqrt(S2[2]))

y.est <- A$pi[1] * dnorm(t, A$Mu[1], sqrt(A$S2[1])) +

A$pi[2] * dnorm(t, A$Mu[2], sqrt(A$S2[2]))

points(t, y, type = "l")

points(t, y.est, type = "l", col = 2, lty = 2)

# assign observations to components - clustering

d <- function(x) which(x == max(x))

apply(A$tau, 1, d)

apply(A$tau, 1, which.max)

# assess misclassification

table(apply(A $tau, 1, which.max), c(rep(1, n1), rep(2, n2)))