涉及 3 个硬币的 EM 问题

Question

我正在使用 EM 算法解决估计问题。问题如下：

您有 3 枚硬币，正面朝上的概率分别为 P1、P2 和 P3。你掷硬币 1。如果硬币 1=H，那么你掷硬币 2；如果硬币 1=T，那么你掷硬币 3。你只记录硬币 2 或 3 是正面还是反面，而不是哪个硬币被翻转。所以观察是一连串的正面和反面，但没有别的。问题是估计 P1、P2 和 P3。

我的 R 代码如下。它不起作用，我不知道为什么。任何想法都将不胜感激，因为我认为这是一个非常狡猾的问题。

本

###############
#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

and of a tail is a minus that, i.e. P1*(1-P2) + (1-P1)*(1-P3) in your example 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

您随后可以做的任何事情都无法区分这些。