Я написал ЭМ для модели смеси, которая состоит из двух олов.
Но это показывает:
"Ошибка: неожиданная '}' в"} "".
Я проверил свой код и не обнаружил проблем с "}".
Буду признателен, если будут обнаружены другие ошибки в моем коде.
Бумага, о которой я говорил, здесь .
И мой код выглядит следующим образом:
#prepare data
slope1 <- -.3;slope2 <- .3;slope3 <- 1.8; slope4 <- 0.5;intercept1 <- 1.5
age <- sample(seq(-2,2,len=201), 40)
grade <- sample(seq(-2,2,len=201), 40)
not_smsa <- sample(seq(-2,2,len=201), 40)
wage <- intercept1 + slope1*age +slope2*grade + slope3*not_smsa + rnorm(length(age),0,.15)
X <- cbind(1, age , grade , not_smsa )
Y <- wage
N<- length(Y)
mydata <- cbind.data.frame(X,Y)
library(mclust)
em <- Mclust(Y)
summary(em, parameters = T)
ans <- lm(wage ~ age + grade + not_smsa ,
data = mydata)
init_params <- c(coef(ans) , coef(ans))
params <-init_params
pi1 <- 0.373321
pi2 <- 0.626679
loglik<- rep(NA, 1000)
mu1 <- X %*% params[1:4]
sigma1 <- sqrt((t((Y-X %*% params[1:4]))%*%(Y-X %*% params[1:4]))/N)
mu2 <- X %*% params[5:8]
sigma2 <- sqrt((t((Y-X %*% params[5:8]))%*%(Y-X %*% params[5:8]))/N)
loglik[1] <- 0
loglik[2] <- sum(log(pi1*dnorm(Y,mu1,sigma1)+pi2*dnorm(Y,mu2,sigma2)))
k<-2
wls <- function(X,Y,W) {
solve(t(X) %*% W %*% X) %*% t(X) %*% W %*% Y
}
while (abs(loglik[k]-loglik[k-1]) >= 0.00001)(na.rm = T){
# E step
w1<-pi1*dnorm(Y,mu1,sigma1)/(pi1*dnorm(Y,mu1,sd=sigma1)+pi2*dnorm(Y,mu2,sigma2))
w2<-pi2*dnorm(Y,mu2,sigma2)/(pi1*dnorm(Y,mu1,sd=sigma1)+pi2*dnorm(Y,mu2,sigma2))
w1[is.na(w1)] <- 0.5
w2[is.na(w2)] <- 0.5
# M step
pi1<-sum(w1)/length(Y)
pi2<-sum(w2)/length(Y)
w <- cbind(w1,w2)
params[1:4] <- wls(X,Y,diag(w[,1]))
params[5:8] <- wls(X,Y,diag(w[,2]))
r1<-Y-X %*% params[1:4]
r2<-Y-X %*% params[5:8]
sigma1h <- sum(w1*r1)/sum(w1)
sigma1h1 <- exp(sigma1h)
sigma2h <- sum(w2*r2)/sum(w2)
sigma2h1 <- exp(sigma2h)
# loglik[k]
loglik[k+1] <- sum(log(pi1*dnorm(Y,mu1,sigma1h1)+pi2*dnorm(Y,mu2,sigma2h1)),na.rm=T)
k<- k+1
}
Спасибо за ваше терпение, и я с нетерпением жду ваших отзывов.