EM涉及3個硬幣的問題

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

和的尾部是一個減去，即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 

你以後可以做什麼都不會區分這些。