# EM algorithm for estimating the mixing proportion of a distribution
# for a pair (x1,x2) of binary observations, when the pairs of observations
# are modelled as being independently drawn from a mixture distribuiton.
# The mixture gives probability 1-p to both x1 and x2 being independently
# from the Bernoulli(0.1) distribution, and probability p of x1 and x2
# being independently from the Bernoulli(0.9) distribution.
#
# The arguments are the data vectors x1 and x2 (binary, equal length), the
# initial value for p, and the number of EM iterations to do.
binmix_em <- function (x1, x2, p, iters)
{
for (i in 1:iters) {
# Find the probabilities for each observation to be generated from
# the two components of the mixture.
p0 <- (1-p) * 0.1^(x1+x2) * 0.9^(2-x1-x2)
p1 <- p * 0.9^(x1+x2) * 0.1^(2-x1-x2)
# Find the conditional probabilities for each observation to come from
# the Bernoulli(0.9) component of the mixture.
r <- p1 / (p0+p1)
# Restimate the mixing proportion parameter.
p <- mean(r)
# Print the log likelihood.
cat(i,p,sum(log(p0+p1)),"\n")
}
# Return the final parameter estimate.
p
}