R语言混合泊松分布EM算法参数振荡收敛问题求助

问题背景 #

用R实现混合泊松分布的EM算法，目标估计参数p、lambda1、lambda2，但迭代中三个参数在两个数值间持续振荡，无法收敛。实现代码及测试数据如下：

原实现代码 #

# E步函数
E_step <- function(Y, lambda1, lambda2, p) {
  poisson1_prob <- dpois(Y, lambda1) * p
  poisson2_prob <- dpois(Y, lambda2) * (1 - p)
  total_poisson_prob <- poisson1_prob + poisson2_prob
  expected_value <- poisson1_prob / total_poisson_prob
  return(expected_value)
}

# M步函数
M_step <- function(X, expected_value) {
  p <- mean(expected_value)
  lambda1 <- sum(expected_value * X) / sum(expected_value)
  lambda2 <- sum((1 - expected_value) * X) / sum(1 - expected_value)
  return(list(p = p, lambda1 = lambda1, lambda2 = lambda2))
}

# 完整EM迭代函数
fullEM <- function(Y, X, lambda1_init, lambda2_init, p_init, num_iter) {
  lambda1 <- lambda1_init
  lambda2 <- lambda2_init
  p <- p_init
  p_new <- numeric(num_iter)
  lambda1_new <- numeric(num_iter)
  lambda2_new <- numeric(num_iter)
  for (i in 1:num_iter) {
    expected_value <- E_step(Y, lambda1, lambda2, p)
    parameters_final <- M_step(X, expected_value)
    p <- parameters_final$p
    lambda1 <- parameters_final$lambda1
    lambda2 <- parameters_final$lambda2
    p_new[i] <- p
    lambda1_new[i] <- lambda1
    lambda2_new[i] <- lambda2
  }
  return(list(p = p_new, lambda1 = lambda1_new, lambda2 = lambda2_new))
}

测试代码及振荡结果 #

# 测试数据
Y <- c(200, 183, 210, 167, 143, 86, 33, 9, 5)
X <- c(0:8)
lambda1_init <- 3
lambda2_init <- 2.3
p_init <- 0.5
num_iter <- 50
result <- fullEM(Y, X, lambda1_init, lambda2_init, p_init, num_iter)
result

运行后参数在两组值间交替：

• p在0.87左右和0.12左右振荡
• lambda1和lambda2在3.44和7.78之间互换

问题原因 #

混合泊松模型存在组分可交换性：模型中两个泊松组分没有固有标签，参数组合(p, λ₁, λ₂)和(1-p, λ₂, λ₁)对应的似然值完全相同，EM算法无法区分哪个组分对应哪个参数，因此会在这两个等价解之间来回跳跃。

解决方案 #

加入参数约束，强制两个组分的lambda值有明确的大小关系（比如λ₁ > λ₂），每次M-step更新后检查约束，若违反则交换lambda值并调整p为1-p，确保参数始终落在同一等价类中，避免振荡。

修改后的代码 #

# E步函数（不变）
E_step <- function(Y, lambda1, lambda2, p) {
  poisson1_prob <- dpois(Y, lambda1) * p
  poisson2_prob <- dpois(Y, lambda2) * (1 - p)
  total_poisson_prob <- poisson1_prob + poisson2_prob
  expected_value <- poisson1_prob / total_poisson_prob
  return(expected_value)
}

# M步函数（不变）
M_step <- function(X, expected_value) {
  p <- mean(expected_value)
  lambda1 <- sum(expected_value * X) / sum(expected_value)
  lambda2 <- sum((1 - expected_value) * X) / sum(1 - expected_value)
  return(list(p = p, lambda1 = lambda1, lambda2 = lambda2))
}

# 完整EM迭代函数（加入参数约束）
fullEM <- function(Y, X, lambda1_init, lambda2_init, p_init, num_iter, tol = 1e-8) {
  lambda1 <- lambda1_init
  lambda2 <- lambda2_init
  p <- p_init
  p_new <- numeric(num_iter)
  lambda1_new <- numeric(num_iter)
  lambda2_new <- numeric(num_iter)
  
  for (i in 1:num_iter) {
    # E步
    expected_value <- E_step(Y, lambda1, lambda2, p)
    # M步
    params <- M_step(X, expected_value)
    # 应用参数约束：确保lambda1 >= lambda2
    if (params$lambda1 < params$lambda2) {
      # 交换lambda并调整p
      new_lambda1 <- params$lambda2
      new_lambda2 <- params$lambda1
      new_p <- 1 - params$p
    } else {
      new_lambda1 <- params$lambda1
      new_lambda2 <- params$lambda2
      new_p <- params$p
    }
    # 检查收敛：参数变化小于阈值则提前停止
    if (abs(new_p - p) < tol && abs(new_lambda1 - lambda1) < tol && abs(new_lambda2 - lambda2) < tol) {
      # 截断存储向量
      p_new <- p_new[1:i]
      lambda1_new <- lambda1_new[1:i]
      lambda2_new <- lambda2_new[1:i]
      break
    }
    # 更新参数
    p <- new_p
    lambda1 <- new_lambda1
    lambda2 <- new_lambda2
    # 存储参数
    p_new[i] <- p
    lambda1_new[i] <- lambda1
    lambda2_new[i] <- lambda2
  }
  return(list(p = p_new, lambda1 = lambda1_new, lambda2 = lambda2_new))
}

测试修改后的代码 #

# 测试数据
Y <- c(200, 183, 210, 167, 143, 86, 33, 9, 5)
X <- c(0:8)
lambda1_init <- 3
lambda2_init <- 2.3
p_init <- 0.5
num_iter <- 50
result <- fullEM(Y, X, lambda1_init, lambda2_init, p_init, num_iter)
result

运行后参数会快速收敛到稳定值，不再振荡：

• p收敛到约0.8719658
• lambda1收敛到约7.779702
• lambda2收敛到约3.445011

额外优化建议 #

• 用对数似然判断收敛：计算每轮的对数似然值，当似然变化小于阈值时停止，比参数变化更可靠。
• 初始化时让lambda1和lambda2差异足够大，减少初始迭代的波动。

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