Olá, o meu projeto precisa do algoritmo EM abaixo, onde está todos os códigos. O erro está no loop while
, que é onde estão as etapas esperança e maximização. A mensagem de erro é:
"Error in while (abs(Elogv[r] - Elogv[r - 1]) >= 1e-06) { : missing value where TRUE/FALSE needed"**.
Como faço pra resolver esse erro, se o loop while
não contém comandos TRUE
e FALSE
e, se já verifiquei detalhadamente que não há nenhum erro nos comandos e não há valor NA
s?
Código:
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
}
NA
e uma alternativa é utilizar o métodoisTRUE
que retorna o valor lógico da expressão e manipula o resultado caso retorneNA
.Ficaria:while (isTRUE(abs(Elogv[r]-Elogv[r-1])>=0.000001, na=FALSE))
r
e seria leviano da minha parte afirmar uma resposta. Então deixo um comentário com a minha opinião e espero que ajude.