在 MATLAB 中测试单峰（单峰）或双峰（双峰）分布

Question

MATLAB 中有没有办法检查直方图分布是单峰还是双峰？

编辑

你认为Hartigan's Dip Statistic 会起作用吗？我尝试将图像传递给它，并获取值0。这是什么意思？

而且，当通过图像时，它是否测试图像的直方图在灰度级别上的分布？

谢谢。

Answer 1

这是一个使用 Nic Price 的 Hartigan's Dip Test 实现来识别单峰分布的脚本。棘手的一点是计算xpdf，这不是概率密度函数，而是排序后的样本。

p_value 是假设原假设为真，获得至少与实际观察到的统计量一样极端的检验统计量的概率。在这种情况下，零假设是分布是单峰的。

close all; clear all;

function [x2, n, b] = compute_xpdf(x)
  x2 = reshape(x, 1, prod(size(x)));
  [n, b] = hist(x2, 40);
  % This is definitely not probability density function
  x2 = sort(x2);
  % downsampling to speed up computations
  x2 = interp1 (1:length(x2), x2, 1:1000:length(x2));
end

nboot = 500;
sample_size = [256 256];

% Unimodal
sample2d = normrnd(0.0, 10.0, sample_size);

[xpdf, n, b] = compute_xpdf(sample2d);
[dip, p_value, xlow, xup] = HartigansDipSignifTest(xpdf, nboot); 

figure;
subplot(1,2,1);
bar(n, b)
title(sprintf('Probability of unimodal %.2f', p_value))

% Bimodal
sample2d = sign(sample2d) .* (abs(sample2d) .^ 0.5);

[xpdf, n, b] = compute_xpdf(sample2d);
[dip, p_value, xlow, xup] = HartigansDipSignifTest(xpdf, nboot); 

subplot(1,2,2);
bar(n, b)
title(sprintf('Probability of unimodal %.2f', p_value))

print -dpng modality.png

Answer 2

有许多不同的方法可以满足您的要求。在最字面意义上，“双峰”意味着有两个峰。但通常，您希望“两个峰”以合理的距离分开，并且您希望它们每个都包含总计数的合理比例。只有您知道什么是您的情况“合理”，但以下方法可能会有所帮助。

1. 创建强度直方图
2. 用cumsum形成累积分布
3. 对于分布（25%、30%、50%，...）之间的不同“切点”值，计算两个分布（切点上方和下方）的均值和标准差。
4. 计算均值之间的距离除以两个分布的标准差之和
5. 该数量将在“最佳切割”时达到最大值

您必须确定该数量的大小对您来说代表“双峰”。这是一些演示我在说什么的代码。它生成不同严重程度的双峰分布 - 两个高斯分布，它们之间的增量增加（步长 = 标准偏差的大小）。我计算了上述数量，并将其绘制为delta 的一系列不同值。然后我在对应于整个分布的 +-1 sigma 的范围内通过这条曲线拟合抛物线。如您所见，当分布变得更加双峰时，会发生两件事：

1. 这条曲线的曲率翻转（从谷到峰）
2. 最大增加（对于高斯大约为 1.33）。

您可以查看您自己的一些分布的这些数量，并决定您希望在哪里设置截止值。

% test for bimodal distribution
close all
for delta = 0:10:50
    a1 = randn(100,100) * 10 + 25;
    a2 = randn(100,100) * 10 + 25 + delta;
    a3 = [a1(:); a2(:)];
    [h hb] = hist(a3, 0:100);
    cs = cumsum(h);
    llimi = find(cs < 0.2 * max(cs(:)));
    ulimi = find(cs > 0.8 * max(cs(:)));
    llim = hb(llimi(end));
    ulim = hb(ulimi(1));
    cuts = linspace(llim, ulim, 20);
    dmean = mean(a3);
    dstd = std(a3);
    for ci = 1:numel(cuts)
        d1 = a3(a3<cuts(ci));
        d2 = a3(a3>=cuts(ci));
        m(ci,1) = mean(d1);
        m(ci, 2) = mean(d2);
        s(ci, 1) = std(d1);
        s(ci, 2) = std(d2);
    end
    q = (m(:, 2) - m(:, 1)) ./ sum(s, 2);
    figure; 
    plot(cuts, q);
    title(sprintf('delta = %d', delta))
    % compute curvature of plot around mean:
    xlims = dmean + [-1 1] * dstd;
    indx = find(cuts < xlims(2) && cuts > xlims(1));
    pf = polyfit(cuts(indx), q(indx), 2);
    m = polyval(pf, dmean);
    fprintf(1, 'coefficients: a = %.2e, peak = %.2f\n', pf(1), m);
end

输出值：

coefficients: a = 1.37e-03, peak = 1.32
coefficients: a = 1.01e-03, peak = 1.34
coefficients: a = 2.85e-04, peak = 1.45
coefficients: a = -5.78e-04, peak = 1.70
coefficients: a = -1.29e-03, peak = 2.08
coefficients: a = -1.58e-03, peak = 2.48

样图：

而 delta = 40 的直方图：