# file "REM3" -- simple EM algorithm for a Poisson model with missing data # data (see Thisted, pp. 188 & 242): count = c(3062, 587, 284, 103, 33, 4, 2) # Here count[i] = the count for i-1 N = sum(count) len = length(count) # initial values: xi = 0.75 lambda = 0.40 cat("i =", 0, " xi =", xi, " lambda =", lambda, "\n") numits = 20 for (i in 1:numits) { # E-step: z = xi / (xi + (1-xi)*exp(-lambda)) * count[1] # M-step: xi = z / N lambda = sum( (0:(len-1)) * count ) / (N-z) # output as we go: cat("i =", i, " z =", z, " xi =", xi, " lambda =", lambda, "\n") }