title: EM Problem involving 3 coins
body: <p>I'm working on an estimation problem using the EM algorithm. The problem is as follows: </p>
  
  <p>You have 3 coins with probabilities of being heads P1, P2, and P3 respectively. You flip coin 1. If coin 1=H, then you flip coin 2; if coin 1=T, then you flip coin 3. You only record whether coin 2 or 3 is heads or tails, not which coin was flipped. So the observations are strings of heads and tails, but nothing else. The problem is to estimate P1, P2, and P3.</p>
  
  <p>My R code to do this is below. It's not working, and I can't figure out why. Any thoughts would be appreciated as I think this is a pretty crafty problem.</p>
  
  <p>Ben</p>
  
  <pre><code>###############
  #simulate data
  p1&lt;-.8
  p2&lt;-.8
  p3&lt;-.3
  tosses&lt;-1000
  rbinom(tosses,size=1,prob=p1)-&gt;coin.1
  pv&lt;-rep(p3,tosses)
  pv[coin.1==1]&lt;-p2
  #face now contains the probabilities of a head
  rbinom(tosses,size=1,prob=pv)-&gt;face
  rm(list=(ls()[ls()!="face"])) 
  #face is all you get to see!
  
  ################
  #e-step
  e.step&lt;-function(x,theta.old) {
      fun&lt;-function(p,theta.old,x) {
          theta.old[1]-&gt;p1
          theta.old[2]-&gt;p2
          theta.old[3]-&gt;p3
          log(p1*p2^x*(1-p2)^(1-x))*(x*p1*p2+(1-x)*p1*(1-p2))-&gt;tmp1 #this is the first part of the expectation
          log((1-p1)*p3^x*(1-p3)^(1-x))*(x*(1-p1)*p3+(1-x)*(1-p1)*(1-p3))-&gt;tmp2 #this is the second
          mean(tmp1+tmp2)
      }
      return(fun)
  }
  #m-step
  m.step&lt;-function(fun,theta.old,face) {
      nlminb(start=runif(3),objective=fun,theta.old=theta.old,x=face,lower=rep(.01,3),upper=rep(.99,3))$par
  }
  
  #initial estimates
  length(face)-&gt;N
  iter&lt;-200
  theta&lt;-matrix(NA,iter,3)
  c(.5,.5,.5)-&gt;theta[1,]
  for (i in 2:iter) {
      e.step(face,theta[i-1,])-&gt;tmp
      m.step(tmp,theta[i-1,],face)-&gt;theta[i,]
      print(c(i,theta[i,]))
      if (max(abs(theta[i,]-theta[i-1,]))&lt;.005) break("conv")
  }
  #note that this thing isn't going anywhere!
  </code></pre>