Matlab / Python中数据的一维高斯混合拟合

Question

我有一条离散曲线y=f(x)。我知道山峰的位置和振幅。我想通过在每个峰值上拟合一个高斯曲线来近似曲线。我应该如何找到优化的高斯参数？我想知道是否有任何内置的函数可以让我的任务变得更简单。

编辑

我固定了高斯的均值，并尝试使用matlab中的lsqcurvefit()对sigma进行优化。MSE较少。然而，我有一个额外的硬约束，即近似曲线的值应该等于峰值处的原始函数。我的模型没有满足这个约束。我在这里粘贴当前的工作代码。我希望有一个解决方案，在峰值遵守硬约束，并在其他点近似拟合曲线。基本思想是近似曲线的参数较少，但仍与原始曲线非常相似。

fun = @(x,xdata)myFun(x,xdata,pks,locs); %pks,locs are the peak locations and amplitudes already available
x0=w(1:6)*0.25; % my initial guess based on domain knowledge 

[sigma resnorm] = lsqcurvefit(fun,x0,xdata,ydata); %xdata and ydata are the original curve data points
recons = myFun(sigma,xdata,pks,locs);
figure;plot(ydata,'r');hold on;plot(recons);

    function f=myFun(sigma,xdata,a,c)
       % a is constant , c is mean of individual gaussians
       f=zeros(size(xdata));
       for i = 1:6 %use 6 gaussians to approximate function
          f = f + a(i) * exp(-(xdata-c(i)).^2 ./ (2*sigma(i)^2));      
       end           
    end

Answer 1

如果你知道你的峰值位置和振幅，那么你所要做的就是找到每个高斯的宽度。你可以把它看作是一个优化问题。

假设你有x和y，它们是你想要近似的曲线的样本。

首先，定义一个函数g()，它将为给定值的宽度构造近似值。g()接受一个参数向量sigma，其中包含每个高斯的宽度。高斯的位置和振幅将被限制在你已经知道的值上。g()输出x中每个点的高斯和近似值。

现在，定义一个以sigma为输入的损失函数L()。L(sigma)返回一个标量来测量误差--给定的近似值(使用sigma)与您试图近似的曲线相差有多大。平方误差是曲线拟合的常见损失函数：

L(sigma) = sum((y - g(sigma)) .^ 2)

现在的任务是搜索sigma的可能值，并找到误差最小的选择。这可以使用各种优化例程来完成。

如果您有Mathworks优化工具箱，则可以使用函数lsqnonlin() (在这种情况下，您不必自己定义L() )。曲线拟合工具箱可能是另一种选择。否则，您可以使用开放源码优化例程(查看cvxopt)。

有几件事需要注意。您需要施加约束，即sigma中的所有值都大于零。您可以告诉优化算法有关此约束的信息。此外，您还需要指定参数的初始猜测(即sigma)。在这种情况下，您可以通过查看每个峰值附近的曲线来选择合理的值。这可能是这样的情况(当损失函数是非凸的)，最终的解决方案是不同的，取决于初始猜测(即，您收敛到局部最小值)。有许多花哨的技术可以处理这种情况，但一件简单的事情就是尝试多个不同的初始猜测，并选择最好的结果。

编辑后添加：

在python中，您可以在scipy.optimize模块中使用优化例程，例如curve_fit()。

编辑2(对编辑的问题的回答)：

如果你的高斯分布彼此有很多重叠，那么求它们的和可能会导致峰值的高度与你的已知值不同。在这种情况下，您可以采用加权总和，并将权重视为另一个要优化的参数。

如果您希望峰值高度恰好等于某些指定值，则可以在优化问题中强制执行此约束。lsqcurvefit()无法做到这一点，因为它只处理参数上的绑定约束。看看fmincon()吧。

Answer 2

您可以使用Expectation–Maximization algorithm在您的数据上拟合混合高斯分布。它不关心数据维度。在MATLAB的文档中，您可以查找gmdistribution.fit或fitgmdist。