如何用Julia拟合一维对角高斯混合模型并约束均值为零？

问题背景 #

你生成了1000个来自两个均值为0、方差分别为1和2的正态混合分布的样本，但使用GaussianMixtures包的默认kmeans初始化拟合GMM后，得到的分量均值不为0——这是因为默认的EM算法会同时优化均值、方差和权重，没有约束均值固定为0的选项。下面提供两种解决方案，确保拟合出的GMM所有分量均值都为0。

方案一：手动实现带均值约束的EM算法 #

对于一维对角GMM，固定均值为0后，只需要优化混合权重和分量方差，EM流程会简化很多：

1. 数据生成（和你的代码一致） #

using Distributions, Random

Random.seed!(42) # 固定种子保证结果可复现
n_samples = 1000
mix_weights = [0.5, 0.5]
components = [Normal(0, 1), Normal(0, sqrt(2))]
mix_model = MixtureModel(components, mix_weights)
samples = rand(mix_model, n_samples)

2. 约束均值的EM实现 #

function fit_zero_mean_gmm(samples, n_components; max_iter=100, tol=1e-6)
    n = length(samples)
    # 初始化：权重均匀分配，方差用样本方差的差异化初始值
    weights = fill(1/n_components, n_components)
    vars = [var(samples) * (i+0.5) for i in 0:n_components-1]
    prev_loglik = -Inf

    for iter in 1:max_iter
        # E步：计算每个样本属于各分量的后验概率（责任度）
        responsibilities = zeros(n, n_components)
        for i in 1:n
            x = samples[i]
            total_density = sum(w * pdf(Normal(0, sqrt(v)), x) for (w, v) in zip(weights, vars))
            for j in 1:n_components
                responsibilities[i,j] = weights[j] * pdf(Normal(0, sqrt(vars[j])), x) / total_density
            end
        end

        # 计算对数似然，判断收敛
        current_loglik = sum(log(sum(weights[j] * pdf(Normal(0, sqrt(vars[j])), x) for j in 1:n_components)) for x in samples)
        if abs(current_loglik - prev_loglik) < tol
            break
        end
        prev_loglik = current_loglik

        # M步：固定均值为0，更新权重和方差
        weights = vec(mean(responsibilities, dims=1))
        for j in 1:n_components
            total_resp = sum(responsibilities[:,j])
            # 均值为0时，方差 = 加权样本平方和 / 总责任度
            vars[j] = sum(responsibilities[i,j] * samples[i]^2 for i in 1:n) / total_resp
        end
    end

    return (weights=weights, variances=vars)
end

# 拟合2分量GMM
result = fit_zero_mean_gmm(samples, 2)
println("混合权重：", result.weights)
println("分量方差：", result.variances)

运行后，你会得到接近真实值（权重[0.5,0.5]，方差[1,2]）的结果，且所有分量均值固定为0。

方案二：修改GaussianMixtures包的EM流程 #

如果你想沿用GaussianMixtures包，可以手动修改EM迭代逻辑，固定均值不更新：

using GaussianMixtures

# 初始化GMM，强制所有分量均值为0
gmm = GMM(2, samples; means=zeros(1, 2)) # 一维数据，2分量，均值初始化为0

function constrained_em!(gmm, samples; max_iter=100, tol=1e-6)
    prev_loglik = -Inf
    n = length(samples)
    for iter in 1:max_iter
        # 复用包内的E步计算责任度
        responsibilities = GaussianMixtures.e_step(gmm, samples)
        current_loglik = GaussianMixtures.loglikelihood(gmm, samples)

        # 收敛判断
        if abs(current_loglik - prev_loglik) < tol
            break
        end
        prev_loglik = current_loglik

        # M步：仅更新权重和方差，不修改均值
        # 更新混合权重
        gmm.w = vec(mean(responsibilities, dims=1))
        # 更新对角方差（一维下为标量）
        for j in 1:gmm.k
            total_resp = sum(responsibilities[:,j])
            cov_val = sum(responsibilities[i,j] * samples[i]^2 for i in 1:n) / total_resp
            gmm.Σ[j] = [cov_val] # 包内协方差存储为矩阵，一维是1x1矩阵
        end
    end
    return gmm
end

# 执行约束EM拟合
constrained_em!(gmm, samples)
println("混合权重：", gmm.w)
println("分量均值：", gmm.μ) # 输出应为[[0.0], [0.0]]
println("分量方差：", [gmm.Σ[j][1,1] for j in 1:2])

这种方法利用了GaussianMixtures的现有工具函数，同时强制均值保持为0，避免了从头实现全部逻辑。

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