如何使用R计算高斯混合模型（Gaussian Mixture Model）的Fisher信息矩阵

计算逻辑 #

首先明确k分量高斯混合模型（GMM）的待估参数集为：各分量权重$\pi_1,\dots,\pi_k$（满足$\sum_{j=1}^k\pi_j=1$）、各分量均值$\mu_1,\dots,\mu_k$、各分量协方差矩阵$\Sigma_1,\dots,\Sigma_k$。
Fisher信息矩阵的计算可按以下步骤进行：

• 先写出单个样本的对数似然：$l(\theta) = \log\left(\sum_{j=1}^k \pi_j \cdot \phi(x;\mu_j,\Sigma_j)\right)$，其中$\phi(\cdot)$为多元高斯分布的概率密度函数，$\theta$为所有参数拼接成的向量
• 求对数似然对$\theta$的一阶偏导，得到得分函数$s(\theta) = \nabla_\theta l(\theta)$
• 最终Fisher信息矩阵为得分函数外积的期望：$I(\theta) = E\left[s(\theta)\cdot s(\theta)^T\right]$，也可通过负的对数似然二阶偏导的期望计算，两种形式等价。

实际应用中，如果是有分量标签的完整数据，可直接推导出解析形式的Fisher矩阵；如果是无标签的观测数据，一般用蒙特卡洛采样近似计算期望，或通过数值求导得到近似结果。

R语言实现方法 #

目前有现成的R包可直接实现该计算，不需要手动推导全量公式：

1. 优先使用mixtools包的gmmpari函数
   该函数可直接返回GMM参数的渐近协方差矩阵，即Fisher信息矩阵的逆，对结果求逆即可得到Fisher信息矩阵，调用示例如下：

# 安装加载包
install.packages("mixtools")
library(mixtools)
library(mvtnorm)

# 模拟示例数据
set.seed(123)
x <- rbind(rmvnorm(500, mean = c(0,0), sigma = diag(2)),
           rmvnorm(500, mean = c(3,3), sigma = diag(2)*1.5))
# 拟合2分量GMM
gmm_fit <- mvnormalmixEM(x, k = 2, verb = FALSE)
# 计算参数渐近协方差（Fisher矩阵的逆）
asy_cov <- gmmpari(gmm_fit)
# 求逆得到Fisher信息矩阵
fim <- solve(asy_cov)

2. 自定义数值求导方案
   如果需要灵活适配不同的参数约束，可搭配numDeriv包的数值求导功能实现，核心逻辑是先定义GMM的对数似然函数，再通过海森矩阵求负得到Fisher矩阵的近似，代码示例如下：

library(numDeriv)
# 定义GMM对数似然函数
gmm_ll <- function(theta, data, k, d) {
  # 拆分参数：权重、均值、协方差
  pi_vec <- c(theta[1:(k-1)], 1 - sum(theta[1:(k-1)]))
  mu_mat <- matrix(theta[k:(k + k*d - 1)], nrow = k, ncol = d)
  sigma_list <- list()
  idx <- k + k*d
  for (j in 1:k) {
    sigma_vec <- theta[idx:(idx + d*(d+1)/2 - 1)]
    sigma_list[[j]] <- matrix(0, nrow = d, ncol = d)
    sigma_list[[j]][upper.tri(sigma_list[[j]], diag = TRUE)] <- sigma_vec
    sigma_list[[j]] <- sigma_list[[j]] + t(sigma_list[[j]]) - diag(diag(sigma_list[[j]]))
    idx <- idx + d*(d+1)/2
  }
  # 计算对数似然
  ll <- sum(log(sapply(1:nrow(data), function(i) {
    sum(sapply(1:k, function(j) pi_vec[j] * dmvnorm(data[i,], mu_mat[j,], sigma_list[[j]])))
  })))
  return(ll)
}
# 代入已得到的参数向量、数据、分量数、样本维度计算海森矩阵
# 示例中theta为你拟合得到的GMM参数拼接后的向量
hess_mat <- hessian(gmm_ll, x = theta, data = x, k = 2, d = 2)
fim <- -hess_mat

内容的提问来源于stack exchange，提问作者Débora