matlab - MATLAB，高斯参数

Question

我必须找到具有至少两个峰值的数据系列的高斯参数。我该如何管理？假设我有yi = f(xi) 并且我需要参数 mu 和 sigma。

我知道我可以取所有数据的对数，然后用 polyfit 计算出来，但这样一来，我就得到了一些我不需要的东西（太长了，说不出原因）。

我应该怎么办？

重要细节：我的 MATLAB 版本没有 normfit。

Answer 1

如果您的 MATLAB 支持kmeans，您可以尝试将数据聚类到两个集群中，然后分别计算每个集群的均值和方差：

%// Cluster bimodal data
idx = kmeans(y, 2);
y1 = y(idx == 1);
y2 = y(idx == 2);

%// Compute means and variances of clusters
M = [mean(y1), mean(y2)];
V = [var(y1), var(y2)];

对于k模式的一般情况，您可以使用以下代码：

idx = kmeans(y, k);    %// Cluster data
C = arrayfun(@(x)y(idx == x), 1:k, 'UniformOutput', false);
M = cellfun(@mean, C); %// Mean of clusters
V = cellfun(@var, C);  %// Variance of clusters

这种方法的好处是它适用于任何数量的集群，只要它是先验已知的。

例子

我们先生成一些任意的双峰高斯数据：

N = 1e4;                    %// Number of samples per mode
M = [1, 5]; V = [0.2, 0.4]; %// Means and variances of two normal distributions
y = bsxfun(@plus, bsxfun(@times, randn(1e4, 1), sqrt(V), M);
y = y(randperm(numel(y)));  %// Shuffle samples

我们应该得到以下直方图：

现在让我们执行 k-means 聚类并计算每个聚类的均值和方差：

idx = kmeans(y, 2);    %// Cluster bimodal data
C = arrayfun(@(x)y(idx == x), 1:k, 'UniformOutput', false);
M = cellfun(@mean, C); %// Mean of clusters
V = cellfun(@var, C);  %// Variance of clusters

我得到的结果是：

M =
    0.9985    4.9802

V =
    0.1949    0.3854

这与原始数据非常接近。

---

如果您没有 MATLAB kmeans，您可以使用 FEX 实现，例如litekmeans.

Answer 2

我已经在这里回答过这种类型的问题几次，每次我都认为“必须可以更简单地做这种事情......”但是，我还没有看到或想到更简单的方法，所以...忍受我:)

如果您事先知道峰值的数量，您可以这样做：

function GaussFit

    % DATA TO REPRODUCE
    mu    = [112  -45];
    sigma = [ 12   24];

    F =[...
        mu(1) + sigma(1)*randn(1e4, 1)
        mu(2) + sigma(2)*randn(1e4, 1)];

    % interpolate with splines through the histogram
    [y,x] = hist(F, 1500);
    G = spline(x,y);

    % Find optimum curve fit 
    P0 = [% mu  S    A
            80  2   2e3;  % (some rough initial estimate)
            -8  12  2e3];        
    P = fminunc(@(P) Obj(P, x,G), P0); % refine the estimate

    % REPRODUCED DATA
    P(:,1:2).'

    figure, clf, hold on
    plot(x, P(1,3)*Gaussian(P(1,1),P(1,2),x) + P(2,3)*Gaussian(P(2,1),P(2,2),x))
    plot(x, ppval(G,x),'r.', 'MarkerSize', 1)

end

% The objective function for the curve fitting optimizer
function val = Obj(P, x,F)

    G = zeros(size(x));    
    for ii = 1:size(P,1);

        mu = P(ii,1);    % mean
        sigma = P(ii,2); % std. deviation
        A = P(ii,3);     % "amplitude"

        G = G + A/sigma/sqrt(2*pi) * exp(-(x-mu).^2/2/sigma^2);

    end    

    val = sum((G-ppval(F,x)).^2);

end

% just a function for plotting
function G = Gaussian(mu,sigma,x)
    G = 1/sigma/sqrt(2*pi) * exp(-(x-mu).^2/2/sigma^2);    
end

结果：

ans =
  112.1633   -45.2013
  12.6777     24.6723

我会说相当不错的结果:)

与往常一样，这种方法有一些缺点。它需要你事先知道

1. 数据集中的峰数
2. 初始估计“足够接近”优化器收敛到实际解决方案

kmeans如果您事先不知道峰值的数量（并且想自动找到峰值的数量） ，则必须使用一些启发式方法来定位数据集中的峰值数量（及其平均值）。

无论如何，重要的是有办法找到峰值的数量，但没有办法自动找到合适的初始估计。如果您只有一个或几十个数据集，仍然可以手动完成初始估计，但除此之外的任何事情都会使上述方法越来越不吸引人。

但是，您可以使用全局优化器，在这种情况下，您不必再提出初始估计。但正是在这一点上，我不禁思考

“这么简单的问题，应该不需要吧！”

但是哦，好吧。