# file "REM2" -- simple EM algorithm for a Poisson model with merged data

# data (simulated from Poisson(2.7):
nstar = 244  # the count for 0 *plus* the count for 1
count = c( 0, 227, 219, 158, 83, 46, 14, 8, 0, 1 )
# Here count[i] = the count for "i", for i >= 2
N = nstar + sum(count)
len = length(count)

# initial value:
lambda = 0.1
cat("i =", 0, "  lambda =", lambda, "\n")

numits = 20

for (i in 1:numits) {

    # E-step:
    z = lambda/(1+lambda) * nstar

    # M-step:
    lambda = (1/N) * ( z + sum((1:len)*count) )

    # output as we go:
    cat("i =", i, "  z =", z, "  lambda =", lambda, "\n")

}