em.norm <- function(Y){
  Yobs <- Y[!is.na(Y)]
  Ymis <- Y[is.na(Y)]
  n <- length(c(Yobs, Ymis))
  r <- length(Yobs)
  # initial values
  mut=1
  sit=0.1
# Define log-likelihood function
  ll <- function(y, mu, sigma2, n){
  -.5*n*log(2*pi*sigma2)-.5*sum((y-mu)^2)/sigma2
}
# Compute the log-likelihood for the initial values,
# and ignoring the missing data mechanism
  lltm1 <- ll(Yobs, mut, sit, n)
  repeat{
  # E-step
  .....
  # M-step
  .....
  # Update parameter values
  .....
  # compute log-likelihood using current estimates,
  # and igoring the missing data mechanism
  llt <- ll(Yobs, mut, sit, n)
  # Print current parameter values and likelihood
  cat(mut, sit, llt, "\n")
  # Stop if converged
  if ( abs(lltm1 - llt) < 0.001) break
  lltm1 <- llt
}
# fill in missing values with new mu.
  return(list(mut=mut,sit=sit))
}
######
em_function.R
   
###
source("em_function.R")
set.seed(123)
x <- rnorm(20,3)
x[16:20] <- NA
x
#[1] 2.439524 2.769823 4.558708 3.070508 3.129288 4.715065 3.460916 1.734939
#[9] 2.313147 2.554338 4.224082 3.359814 3.400771 3.110683 2.444159       NA
#[17]      NA       NA       NA       NA
em.norm(x)
#2.614288 1.393861 -26.84626
#3.01786 0.9029745 -23.04815
#3.118753 0.7293558 -22.09283
#3.143977 0.6827701 -21.88989
#3.150282 0.6709248 -21.8435
#3.151859 0.6679511 -21.83223
#3.152253 0.6672069 -21.82943
#3.152351 0.6670208 -21.82873
#$mut    [1] 3.152351
#$sit      [1] 0.6670208
mean(x,na.rm=T)                                             #[1] 3.152384
mean(x^2,na.rm=T) - mean(x,na.rm=T)^2     # variance [1] 0.6669587

xiexie

有木有二元变量的程序啊，对于一元变量实际上得到的就是观测值的最大似然估计，是吗

很厉害！