如何解决R中循环EM算法中的错误

Question

我的项目需要下面的EM算法，所有的代码在哪里。错误出现在 while 循环 中，这是希望和最大化步骤所在的位置。错误消息是“错误 in while (abs (Elogv [r] - Elogv [r - 1])> = 1e-06) {: 需要 TRUE / FALSE 的缺失值”。如果 while 循环不包含 true 和 false 命令，并且我已经详细检查了命令中没有错误并且没有 NA 的值，我该如何解决这个错误？感谢关注，谁能救救我。

n=100
u<-runif(n)
QUANTIL <- function(u){
  Q <- rep(NA, length(u))
  for (i in 1:length(u)) {
    if(u[i] <  0.2634253829){
      Q[i] <- 1*tan(pi*(0.9490353482*u[i]-0.5))+0
    }
    if(u[i]>=0.2634253829 && u[i] < 0.7365746171){
      Q[i] <-  1*qnorm(1.4428629504*u[i]-0.2214315)+0
    }
    if(u[i]>0.7365746171){
      Q[i] <- 1*tan(pi*(0.9490353482*u[i]-0.4490353))+0
    } 
  }
  return(Q)
}
x<-QUANTIL(u)
y<-c(sort(x))
i<-seq(1,n)
v<-c(i/(n+1))

t<-QUANTIL(v)
mi<-median(y)
s<-c(y[26:73])
sigma<-sqrt(sum((s-mi)^2)/(n-1))
p=0.4731492342

alpha<-(2*t^3)/(1+t^2)^2
beta<-(1-t^2)/(1+t^2)^2
eta<-(t^4-t^2)/(1+t^2)^2
lambda<-2*t/(1+t^2)^2
gama<-(-t^2)
delta<-2*t

k<-((p*0.6930665173/sigma*sqrt(2*pi))*exp((-1/2*sigma^2)*((y-mi)^2)))/(((p*0.6930665173/sigma*sqrt(2*pi))*exp((-1/2*sigma^2)*(y-mi)^2))+((((1-p)*1.0537015317/sigma*pi))*(1/(1+((y-mi)/sigma)^2))))
r<-2
Elogv<-sum(k*((-1/2)*((y-mi)/sigma)^2))-sum(k*log(sigma*sqrt(2*pi)))-sum((1-k)*log(sigma*pi))-sum((1-k)*log(1+((y-mi)/sigma)^2))+sum(k*log(p))+(n-sum(k))*log(1-p)+log(0.6930665173)*sum(k)+log(1.0537015317)*sum(1-k)
Elogv[1]<-0

while (abs(Elogv[r]-Elogv[r-1])>=0.000001) {

  w<-(2*beta-2*k*beta+k)
  q<-k*delta+2*lambda*(1-k)
  sigma<-(sum(y*w)*sum(q)-sum(w)*sum(y*q))/(-2*sum(alpha*(1-k))*sum(q)+sum(w)*sum(k*gama-1)+2*sum(w)*sum(eta*(1-k)))                                                  
  mi<-(sum(y*w)+2*sigma*sum(alpha*(1-k)))/sum(w)
  k<-((p*0.6930665173/sigma*sqrt(2*pi))*exp((-1/2*sigma^2)*((y-mi)^2)))/(((p*0.6930665173/sigma*sqrt(2*pi))*exp((-1/2*sigma^2)*(y-mi)^2))+((((1-p)*1.0537015317/sigma*pi))*(1/(1+((y-mi)/sigma)^2))))
  Elogv[r]<-sum(k*((-1/2)*((y-mi)/sigma)^2))-sum(k*log(sigma*sqrt(2*pi)))-sum((1-k)*log(sigma*pi))-sum((1-k)*log(1+((y-mi)/sigma)^2))+sum(k*log(p))+(n-sum(k))*log(1-p)+log(0.6930665173)*sum(k)+log(1.0537015317)*sum(1-k)
  r<-r+1

Answer 1

在我看来，Elogv 的长度是 1？因此 Elogv[r] 没有条目（r 为 2！），即计算结果为 NA，因此 abs(Elogv[r]-Elogv[r-1]) 为 NA。

在开始循环之前你需要 Elogv[2]