三正态分布混合模型EM算法运行报错，请求技术协助

问题描述 #

需要实现针对三正态分布混合模型的EM算法（均值和方差均未知），数据为500行的Height列（记为S）。编写了负对数似然函数和EM迭代代码，但运行后出现大量错误，需排查解决。

原负对数似然函数代码：

S <- df$Height

neg.logl = function(theta, p1, p2,x){
  mu1 = theta[1]
  mu2 = theta[2]
  mu3 = theta[3]
  sigma1 = exp(theta[4]) 
  sigma2 = exp(theta[5]) 
  sigma3 = exp(theta[6])
  f.x = p1 * dnorm(x, mu1, sigma1) + 
    p2*dnorm(x, mu2, sigma2) +
    (1-p1-p2)*dnorm(x, mu3, sigma3)
  sum(-log(f.x))
}

原EM迭代代码：

set.seed(30000)
max.iter = 10; n=nrow(df)
p1.comp1 = matrix(NA, max.iter, n)
p2.comp1 = matrix(NA, max.iter, n)
theta_EM = matrix(NA, max.iter, 6)
p1.comp1[1, ] = runif(n,0.1,0.3)
p2.comp1[1, ] = runif(n,0.2,0.3)

# Running first M step
theta_EM[1, ] = optim(c(-0.5, 0.5, -0.2, 0.2), function(theta)
  neg.logl(theta, p1.comp1[1,],p2.comp1[1, ], S))$par

# Running alternate E and M steps
for(i in 2:max.iter){
# E step: update probability memberships
p1.temp = cbind(dnorm(S, theta_EM[i-1, 1], exp(theta_EM[i-1, 4])),
  dnorm(S, theta_EM[i-1, 2], exp(theta_EM[i-1, 5])),
  dnorm(S, theta_EM[i-1, 3], exp(theta_EM[i-1, 6])))
p2.temp = cbind(dnorm(S, theta_EM[i-1, 1], exp(theta_EM[i-1, 4])),
  dnorm(S, theta_EM[i-1, 2], exp(theta_EM[i-1, 5])),
  dnorm(S, theta_EM[i-1, 3], exp(theta_EM[i-1, 6])))
p1.comp1[i, ] = p1.temp[ , 1] / rowSums(p1.temp)
p2.comp1[i, ] = p2.temp[ , 1] / rowSums(p2.temp)

# M step: update parameter estimates
theta_EM[i, ] = optim(theta_EM[i-1, ],
function(theta) neg.logl(theta, p1.comp1[i,],p2.comp1[i, ], S))$par
}

#Final
theta.last = theta_EM[max.iter, ]
p1.temp = cbind(dnorm(S, theta.last[1], exp(theta.last[4])),
  dnorm(x, theta.last[2], exp(theta.last[5])))
p2.temp = cbind(dnorm(S, theta.last[2], exp(theta.last[5])),
  dnorm(x, theta.last[2], exp(theta.last[6])))
p3.temp = cbind(dnorm(S, theta.last[3], exp(theta.last[6])),
  dnorm(x, theta.last[2], exp(theta.last[6])))
p1.last = p1.temp[ , 1]/rowSums(p1.temp)
p2.last = p2.temp[ , 1]/rowSums(p2.temp)
p2.last = p3.temp[ , 1]/rowSums(p3.temp)
logl.last = - neg.logl(theta.last, p1.last, p2.last,S)
cat("EP:",
c(theta.last[1], theta.last[2], theta.last[3],exp(theta.last[4]), 
    exp(theta.last[5],exp(theta.last[6])),
    "
logl:", logl.last)

错误排查与修正点 #

1. 初始参数维度不匹配 #

optim初始值向量长度需与theta的6个参数（mu1,mu2,mu3,sigma1_log,sigma2_log,sigma3_log）一致，但原代码仅传入4个值，导致访问theta[5]和theta[6]时出错。

修正：改为长度为6的初始向量，例如c(-0.5, 0.5, 0, 0, 0, 0)。

2. E步骤成员概率计算逻辑错误 #

原代码中p2.temp完全复制p1.temp，且p2.comp1[i,]错误取第一个成分的密度，违背了三成分后验概率的计算逻辑：每个样本的后验概率是该样本在对应成分的密度乘以成分权重后归一化，需同时维护三个成分的后验概率。

修正：用单个矩阵存储三个成分的后验概率，避免拆分存储导致的逻辑混乱。

3. M步骤负对数似然参数传递错误 #

原neg.logl函数的p1,p2应为成分先验权重（标量），但代码中传递的是每个样本的后验概率（向量），完全不符合似然函数的定义。此外，三正态混合的M步骤有解析解，无需用optim数值优化，效率和稳定性更优。

修正：改用M步骤的解析解更新参数，或重新定义负对数似然函数使其接受完整的模型参数。

4. 最终计算部分的变量与语法错误 #

• 未定义变量x，应替换为S；
• p1.temp/p2.temp/p3.temp列数错误，需包含三个成分的密度；
• 重复赋值p2.last覆盖结果，导致第三个成分概率丢失；
• cat语句中exp(theta.last[5],exp(theta.last[6]))是错误的函数调用，且字符串拼接格式错误。

5. 初始后验概率合理性问题 #

原代码随机生成的p1.comp1和p2.comp1未保证p1.comp1 + p2.comp1 <=1，会导致1-p1-p2为负，进而引发密度计算错误。

修正后的完整代码（推荐解析解版本） #

set.seed(30000)
S <- df$Height
n <- length(S)
max.iter <- 10

# 初始化参数：3个成分的均值、标准差、权重（权重和为1）
init_mu <- c(rnorm(3, mean(S), sd(S)))
init_sigma <- rep(sd(S), 3)
init_p <- c(0.3, 0.3, 0.4)

# 存储迭代过程参数
params <- data.frame(iter = 1:max.iter,
                     mu1 = NA, mu2 = NA, mu3 = NA,
                     sigma1 = NA, sigma2 = NA, sigma3 = NA,
                     p1 = NA, p2 = NA, p3 = NA,
                     log_likelihood = NA)
params[1, c("mu1","mu2","mu3")] <- init_mu
params[1, c("sigma1","sigma2","sigma3")] <- init_sigma
params[1, c("p1","p2","p3")] <- init_p

# 对数似然函数
log_likelihood <- function(mu, sigma, p, x) {
  sum(log(p[1]*dnorm(x, mu[1], sigma[1]) + 
            p[2]*dnorm(x, mu[2], sigma[2]) + 
            p[3]*dnorm(x, mu[3], sigma[3])))
}

params$log_likelihood[1] <- log_likelihood(init_mu, init_sigma, init_p, S)

# EM迭代
for (i in 2:max.iter) {
  # E步骤：计算后验概率
  prev_mu <- params[i-1, c("mu1","mu2","mu3")]
  prev_sigma <- params[i-1, c("sigma1","sigma2","sigma3")]
  prev_p <- params[i-1, c("p1","p2","p3")]
  
  dens1 <- prev_p[1] * dnorm(S, prev_mu[1], prev_sigma[1])
  dens2 <- prev_p[2] * dnorm(S, prev_mu[2], prev_sigma[2])
  dens3 <- prev_p[3] * dnorm(S, prev_mu[3], prev_sigma[3])
  total_dens <- dens1 + dens2 + dens3
  
  gamma1 <- dens1 / total_dens
  gamma2 <- dens2 / total_dens
  gamma3 <- dens3 / total_dens
  
  # M步骤：解析解更新参数
  new_mu1 <- sum(gamma1 * S) / sum(gamma1)
  new_mu2 <- sum(gamma2 * S) / sum(gamma2)
  new_mu3 <- sum(gamma3 * S) / sum(gamma3)
  
  new_sigma1 <- sqrt(sum(gamma1 * (S - new_mu1)^2) / sum(gamma1))
  new_sigma2 <- sqrt(sum(gamma2 * (S - new_mu2)^2) / sum(gamma2))
  new_sigma3 <- sqrt(sum(gamma3 * (S - new_mu3)^2) / sum(gamma3))
  
  new_p1 <- mean(gamma1)
  new_p2 <- mean(gamma2)
  new_p3 <- mean(gamma3)
  
  # 存储参数与对数似然
  params[i, c("mu1","mu2","mu3")] <- c(new_mu1, new_mu2, new_mu3)
  params[i, c("sigma1","sigma2","sigma3")] <- c(new_sigma1, new_sigma2, new_sigma3)
  params[i, c("p1","p2","p3")] <- c(new_p1, new_p2, new_p3)
  params[i, "log_likelihood"] <- log_likelihood(c(new_mu1,new_mu2,new_mu3), 
                                                c(new_sigma1,new_sigma2,new_sigma3), 
                                                c(new_p1,new_p2,new_p3), S)
}

# 输出最终结果
cat("最终估计参数：\n")
cat("均值：", tail(params[,c("mu1","mu2","mu3")],1), "\n")
cat("标准差：", tail(params[,c("sigma1","sigma2","sigma3")],1), "\n")
cat("成分权重：", tail(params[,c("p1","p2","p3")],1), "\n")
cat("对数似然：", tail(params$log_likelihood,1), "\n")

# 查看迭代过程
print(params)

关键说明 #

1. 优先使用M步骤的解析解，避免数值优化带来的不稳定和效率问题；
2. 初始参数建议基于数据统计量（如均值、分位数）设置，降低收敛到局部最优的概率；
3. 可添加收敛判断（如对数似然变化小于阈值时停止），替代固定迭代次数。

内容的提问来源于stack exchange，提问作者user11607046