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

#algorithm #r

#алгоритм #r

Вопрос:

Я работаю над проблемой оценки с использованием алгоритма 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!

1. Что именно означает «не работает»? И вам, вероятно, следует добавить тег «домашнее задание» — оценка вероятности подбрасывания монет не имеет большого смысла вне этого (если только вы не бюджетное казино Вегаса ;), то есть …)

2. Я добавил этот тег «домашнее задание». Оценки просто не идут туда, куда они должны идти. Код выполняется, он просто выдает мусор.

Ответ №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

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