MATLAB evalclusters函数的Julia等效实现：Kmeans最优K求解

Julia的生态中没有和MATLAB evalclusters 完全对应的内置函数，但可以通过手动实现Gap统计量结合现有聚类包完成相同功能。以下是具体方案：

核心思路 #

Gap统计量的核心逻辑是：对比真实数据的聚类内离散度（WCSS）与参考分布（如数据边界内的均匀分布）的聚类内离散度，通过两者的对数差值确定最优K。我们将基于Clustering.jl（Julia主流聚类包）实现这一逻辑。

步骤1：安装依赖包 #

先安装所需的Julia包：

using Pkg
Pkg.add(["Clustering", "Distributions", "StatsBase"])

步骤2：实现Gap统计量计算函数 #

以下代码完整实现Gap统计量的计算、参考分布生成及最优K选择：

using Clustering, Distributions, StatsBase

# 计算聚类内加权平方和的对数
function compute_log_wcss(data, assignments, centers)
    wcss = 0.0
    for k in 1:length(centers)
        cluster_points = data[:, assignments .== k]
        if !isempty(cluster_points)
            wcss += sum(colwise(SqEuclidean(), cluster_points, centers[:,k]))
        end
    end
    return log(wcss)
end

# 生成参考分布数据（均匀分布在数据的边界框内）
function generate_reference_data(data, n_samples)
    n_features, n_obs = size(data)
    ref_data = similar(data, n_features, n_samples)
    for i in 1:n_features
        min_val = minimum(data[i,:])
        max_val = maximum(data[i,:])
        ref_data[i,:] = rand(Uniform(min_val, max_val), n_samples)
    end
    return ref_data
end

# 计算Gap统计量
function gap_statistic(data, K_list; B=10)
    n_obs = size(data, 2)
    gaps = zeros(length(K_list))
    gap_stds = zeros(length(K_list))
    
    for (idx, K) in enumerate(K_list)
        # 计算真实数据的log(WCSS)
        kmeans_res = kmeans(data, K; maxiter=1000, init=:kmpp)
        log_wcss_actual = compute_log_wcss(data, kmeans_res.assignments, kmeans_res.centers)
        
        # 多次生成参考分布并计算log(WCSS)
        log_wcss_ref = zeros(B)
        for b in 1:B
            ref_data = generate_reference_data(data, n_obs)
            ref_kmeans = kmeans(ref_data, K; maxiter=1000, init=:kmpp)
            log_wcss_ref[b] = compute_log_wcss(ref_data, ref_kmeans.assignments, ref_kmeans.centers)
        end
        
        # 计算Gap值和标准差
        gaps[idx] = mean(log_wcss_ref) - log_wcss_actual
        gap_stds[idx] = std(log_wcss_ref) * sqrt(1 + 1/B)
    end
    
    return gaps, gap_stds, K_list
end

# 根据Gap统计量规则选择最优K
function find_optimal_K(gaps, gap_stds, K_list)
    # 遵循Tibshirani的规则：找到最小的K，使得gap(K) ≥ gap(K+1) - gap_stds(K+1)
    for i in 1:length(K_list)-1
        if gaps[i] >= gaps[i+1] - gap_stds[i+1]
            return K_list[i]
        end
    end
    # 若未满足条件，返回最大候选K
    return K_list[end]
end

步骤3：使用示例（对应MATLAB原代码） #

假设你的W0是MATLAB格式的样本×特征矩阵，在Julia中需要转置为特征×样本（Clustering.jl的输入要求）：

# 转换数据维度（MATLAB样本×特征 → Julia特征×样本）
data = permutedims(W0)
max_K = 10
K_list = 1:max_K

# 计算Gap统计量
gaps, gap_stds, Ks = gap_statistic(data, K_list)
optimal_K = find_optimal_K(gaps, gap_stds, Ks)

# 用最优K执行K-means得到聚类结果
final_kmeans = kmeans(data, optimal_K; maxiter=1000, init=:kmpp)
idx = final_kmeans.assignments  # 对应MATLAB的evaluation.OptimalY

关键细节说明 #

• 参考分布：代码中生成的参考分布与MATLAB evalclusters 默认逻辑一致（均匀分布在数据的边界框内）。
• 初始化方法：使用kmpp（K-means++）初始化，与MATLAB K-means默认行为对齐。
• 参数调整：B是参考分布的生成次数，默认10次，可根据精度需求增大（但会增加计算时间）。

内容的提问来源于stack exchange，提问作者arrow-economist