ЭМ Проблема с участием 3 монет

Question

Я работаю над проблемой оценки, используя алгоритм EM. Проблема заключается в следующем:

У вас есть 3 монеты с вероятностью быть головами P1, P2 и P3 соответственно. Вы подбрасываете монету 1. Если монета 1 = H, тогда вы подбрасываете монету 2; если монета 1 = T, то вы подбрасываете монету 3. Вы записываете, является ли монета 2 или 3 головами или хвостами, а не какая монета была подброшена. Таким образом, наблюдения представляют собой цепочки голов и хвостов, но не более того. Проблема состоит в том, чтобы оценить P1, P2 и P3.

Мой код R для этого ниже. Это не работает, и я не могу понять, почему. Будем благодарны за любые мысли, так как я думаю, что это довольно коварная проблема.

Ben

###############
#simulate data
p1<-.8
p2<-.8
p3<-.3
tosses<-1000
rbinom(tosses,size=1,prob=p1)->coin.1
pv<-rep(p3,tosses)
pv[coin.1==1]<-p2
#face now contains the probabilities of a head
rbinom(tosses,size=1,prob=pv)->face
rm(list=(ls()[ls()!="face"])) 
#face is all you get to see!

################
#e-step
e.step<-function(x,theta.old) {
    fun<-function(p,theta.old,x) {
        theta.old[1]->p1
        theta.old[2]->p2
        theta.old[3]->p3
        log(p1*p2^x*(1-p2)^(1-x))*(x*p1*p2+(1-x)*p1*(1-p2))->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))->tmp2 #this is the second
        mean(tmp1+tmp2)
    }
    return(fun)
}
#m-step
m.step<-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)->N
iter<-200
theta<-matrix(NA,iter,3)
c(.5,.5,.5)->theta[1,]
for (i in 2:iter) {
    e.step(face,theta[i-1,])->tmp
    m.step(tmp,theta[i-1,],face)->theta[i,]
    print(c(i,theta[i,]))
    if (max(abs(theta[i,]-theta[i-1,]))<.005) break("conv")
}
#note that this thing isn't going anywhere!

Answer 1

Вы не можете оценить P1, P2 и P3 отдельно. Единственная полезная информация - это доля записанных головок и общее количество наборов переворотов (каждый набор переворотов независим, поэтому порядок не имеет значения). Это все равно что пытаться решить одно уравнение для трех неизвестных, а это невозможно сделать.

Вероятность записи головы равна P1*P2 + (1-P1)*P3, что в вашем примере равно 0,7

и хвоста - это один минус, то есть P1*(1-P2) + (1-P1)*(1-P3) в вашем примере 0.3

Вот простой симулятор

#simulate data
sim <- function(tosses, p1, p2, p3) { 
          coin.1   <- rbinom(tosses, size=1, prob=p1)
          coin.2   <- rbinom(tosses, size=1, prob=p2)
          coin.3   <- rbinom(tosses, size=1, prob=p3)
          ifelse(coin.1 == 1, coin.2, coin.3) # returned 
    }

Ниже приведены иллюстрации, все из которых дают 0,7 (с некоторыми случайными колебаниями)

> mean(sim(100000, 0.8, 0.8, 0.3))
[1] 0.70172
> mean(sim(100000, 0.2, 0.3, 0.8))
[1] 0.69864
> mean(sim(100000, 0.5, 1.0, 0.4))
[1] 0.69795
> mean(sim(100000, 0.3, 0.7, 0.7))
[1] 0.69892
> mean(sim(100000, 0.5, 0.5, 0.9))
[1] 0.70054
> mean(sim(100000, 0.6, 0.9, 0.4))
[1] 0.70201

Ничто из того, что вы можете сделать впоследствии, не сможет их различить.