# STA 2102, Fall 200, Assignment #4. # FIND ML ESTIMATES FOR A MIXTURE OF EXPONENTIALS USING EM. The arguments # are the vector of data points, the initial values for the means of the # two components, the initial value for the probability of a data point # coming from the first component, the number of iterations of EM to do, # and a debug option. If the debug option is zero, no debug information # is printed; if it is greater than zero, debug information is printed at # intervals given by its value. The value returned is the vector of parameter # estimates, for the mean time to failure of defective machines, the mean time # to failure of good machines, and the probability a machine is defective. em <- function (x,m1=mean(x)*0.9,m2=mean(x)*1.1,p=0.5,n=20,debug=0) { if (debug!=0) { cat(0,m1,m2,p,log.likelihood(x,m1,m2,p),"\n") } for (i in 1:n) { odds <- (p * dexp(x,1/m1)) / ((1-p) * dexp(x,1/m2)) r <- odds / (1+odds) p <- mean(r) m1 <- sum(r*x) / sum(r) m2 <- sum((1-r)*x) / sum(1-r) if (debug!=0 && i%%debug==0) { cat(i,m1,m2,p,log.likelihood(x,m1,m2,p),"\n") } } return (m1,m2,p) } # LOG LIKELIHOOD FUNCTION FOR THE EXPONENTIAL MIXTURE. Not needed to do EM, # but useful for debugging. log.likelihood <- function (x,m1,m2,p) { sum (log (p*dexp(x,1/m1) + (1-p)*dexp(x,1/m2))) } # FUNCTION TO GENERATE TEST DATA FROM AN EXPONENTIAL MIXTURE. gen <- function (m1,m2,p,N) { w <- runif(N)