在R中拟合截断正态分布混合模型以检测变量模态数

我需要在R中检测1-7量表调查数据的模态数，因此要拟合截断正态分布混合模型，但目前没找到合适的实现方法。

模拟三组分截断正态混合数据的代码 #

# install.packages("truncnorm")
library(truncnorm)

set.seed(123)

# 模拟第一个截断正态分布
mu1 <- 0.5; sd1 <- 0.5
a1 <- 0; b1 <- 7 # 上下界
N1 <- 100
sample1 <- rtruncnorm(N1, a=a1, b=b1, mean=mu1, sd=sd1)

# 模拟第二个截断正态分布
mu2 <- 5; sd2 <- 1.5
a2 <- 0; b2 <- 7
N2 <- 110
sample2 <- rtruncnorm(N2, a=a2, b=b2, mean=mu2, sd=sd2)

# 模拟第三个截断正态分布
mu3 <- 6.8; sd3 <- 0.5
a3 <- 0; b3 <- 7
N3 <- 60
sample3 <- rtruncnorm(N3, a=a3, b=b3, mean=mu3, sd=sd3)

# 合并样本
combined_sample <- c(sample1, sample2, sample3)
data <- data.frame(combined_sample, 
                   class = c(rep("sample1", N1), rep("sample2", N2), rep("sample3", N3))) # 真实类别

# 绘制合并分布
hist(data$combined_sample, breaks=30, main="合并的截断正态分布", 
     xlab="取值", ylab="频数", col="lightblue", border="black")
abline(v=c(mu1, mu2, mu3), col=c("red", "blue", "green"), lwd=2) # 添加真实均值作为参考

无截断高斯混合模型的局限性 #

我常用mclust包拟合无截断的高斯混合模型，但它无法可靠识别截断数据的模态数，会高估组分数量：

# install.packages("mclust")
library(mclust)

# 最低BIC对应的模型有4个组分
mclustBIC(data$combined_sample)

需求 #

请问在R中是否有方法可以适配潜在正态分布的截断特性？

在R中可以通过以下几种方式拟合截断正态分布混合模型，适配数据的截断特性：

方法1：使用flexmix包自定义截断正态组分 #

flexmix支持自定义分布的混合模型，我们可以基于截断正态的密度函数定义组分，然后拟合模型并选择最优组分数量。

# 安装并加载包
# install.packages(c("flexmix", "truncnorm"))
library(flexmix)
library(truncnorm)

# 定义截断正态的密度函数
dtruncnorm_mix <- function(x, mu, sigma, a=0, b=7) {
  dtruncnorm(x, a=a, b=b, mean=mu, sd=sigma)
}

# 定义flexmix的组件
setMethod("FLXMRdensity", signature(model="truncnorm"), 
          function(model, y, x, w, ...) {
            pars <- FLXgetPars(model)
            mu <- pars[,1]
            sigma <- pars[,2]
            logLik <- sum(w * log(dtruncnorm_mix(y, mu=mu, sigma=sigma)))
            new("FLXMRdensity", 
                parameters=pars, 
                logLik=logLik)
          })

# 初始化参数（尝试2-5个组分）
init_pars <- function(k) {
  # 用k-means初始划分均值
  km <- kmeans(data$combined_sample, centers=k)
  mu_init <- km$centers[,1]
  sigma_init <- rep(sd(data$combined_sample), k)
  cbind(mu_init, sigma_init)
}

# 拟合不同组分数量的模型
bic_values <- c()
for(k in 2:5) {
  model <- flexmix(combined_sample ~ 1, data=data, 
                   k=k, 
                   model=FLXMRdensity("truncnorm"),
                   init=list(parameters=init_pars(k)))
  bic_values[k-1] <- BIC(model)
}

# 查看最优组分数量（BIC最小的k）
which.min(bic_values) + 1 # 应返回3，对应真实组分数量

# 获取最优模型的参数
best_k <- which.min(bic_values) + 1
best_model <- flexmix(combined_sample ~ 1, data=data, 
                      k=best_k, 
                      model=FLXMRdensity("truncnorm"),
                      init=list(parameters=init_pars(best_k)))
parameters(best_model)

方法2：使用mixtools包的截断正态混合模型 #

mixtools提供了truncnormmixEM函数，专门用于拟合截断正态分布的混合模型，同时支持选择组分数量。

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

# 拟合截断正态混合模型（指定上下界0-7，尝试3个组分）
set.seed(123)
fit <- truncnormmixEM(data$combined_sample, 
                      lambda=NULL, mu=NULL, sigma=NULL,
                      lower=0, upper=7, k=3)

# 查看拟合结果
fit$lambda # 组分权重
fit$mu     # 组分均值
fit$sigma  # 组分标准差

# 尝试不同k值，通过BIC选择最优模型
bic_trunc <- function(k) {
  fit <- truncnormmixEM(data$combined_sample, lower=0, upper=7, k=k)
  # 计算BIC：-2*logLik + k*(3)*log(n) （每个组分有lambda、mu、sigma三个参数）
  n <- length(data$combined_sample)
  -2*fit$loglik + k*3*log(n)
}

# 计算2-5个组分的BIC
bic_k <- sapply(2:5, bic_trunc)
best_k <- which.min(bic_k) + 1 # 应得到3

方法3：自定义似然函数用optim拟合 #

如果需要更灵活的控制，可以自定义混合模型的对数似然函数，用optim进行参数估计。

# 定义对数似然函数
log_lik_trunc_mix <- function(pars, x, lower=0, upper=7) {
  k <- length(pars) %/% 3
  lambda <- pars[1:k]
  mu <- pars[(k+1):(2*k)]
  sigma <- pars[(2*k+1):(3*k)]
  
  # 确保权重和为1，标准差为正
  lambda <- lambda / sum(lambda)
  sigma <- abs(sigma)
  
  # 计算每个观测的对数似然
  ll <- sapply(x, function(val) {
    sum(lambda * dtruncnorm(val, a=lower, b=upper, mean=mu, sd=sigma))
  })
  -sum(log(ll)) # 负对数似然用于最小化
}

# 初始化参数（k=3）
init_lambda <- c(1/3, 1/3, 1/3)
init_mu <- c(1, 5, 7)
init_sigma <- c(1, 1, 1)
init_pars <- c(init_lambda, init_mu, init_sigma)

# 优化拟合
set.seed(123)
fit_optim <- optim(init_pars, log_lik_trunc_mix, x=data$combined_sample,
                   method="L-BFGS-B", 
                   lower=c(rep(0.01,3), rep(0,3), rep(0.01,3)),
                   upper=c(rep(0.99,3), rep(7,3), rep(5,3)))

# 提取拟合参数
k <- 3
est_lambda <- fit_optim$par[1:k] / sum(fit_optim$par[1:k])
est_mu <- fit_optim$par[(k+1):(2*k)]
est_sigma <- abs(fit_optim$par[(2*k+1):(3*k)])

# 查看结果
est_lambda
est_mu
est_sigma

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