EM算法--隐含变量

# 创建一个示例数据集，假设数据由两个高斯分布混合而成
set.seed(123)
data <- c(rnorm(1000, mean = 2, sd = 1), rnorm(1000, mean = 7, sd = 2))
#概率密度函数
f <- function(data,mu,sigma) {
  f<-1/(sqrt(2*pi)*sigma)*exp(-(data-mu)^2/2*sigma^2)
  return(f)
}
# 定义EM算法函数
em_algorithm <- function(data, k = 2, max_iterations = 100) {
  # 初始化参数
  n <- length(data)
  w <- rep(1/k, k)  # 混合系数
  mu <- sample(data, k)  # 均值参数
  sigma <- rep(1, k)  # 标准差参数
  # EM算法迭代
  for (iteration in 1:max_iterations) {
    # E步骤：计算每个数据点属于每个分布的后验概率
    posterior <- matrix(0, nrow = n, ncol = k)
    for (i in 1:n) {
      for (j in 1:k) {
        posterior[i, j] <- w[j] * f(data[i],mu[j],sigma[j])
      }
      posterior[i, ] <- posterior[i, ] / sum(posterior[i, ])  # 归一化
    }
    
    # M步骤：更新参数估计
    for (j in 1:k) {
      w[j] <- sum(posterior[, j]) / n  # 更新混合系数
      mu[j] <- sum(posterior[, j] * data) / sum(posterior[, j])  # 更新均值参数
      sigma[j] <- sqrt(sum(posterior[, j] * (data - mu[j])^2) / sum(posterior[, j]))  # 更新标准差参数
    }
  }
  # 返回最终的估计结果
  return(list(weights = w, means = mu, standard_deviations = sigma))
}
# 使用EM算法估计参数
result <- em_algorithm(data, k = 2, max_iterations = 100)
result$means
result$standard_deviations