【问题标题】:How to fit a gaussian to data in matlab/octave?如何将高斯拟合到 matlab/octave 中的数据?
【发布时间】:2012-10-28 17:39:59
【问题描述】:

我有一组带有峰值的频率数据,我需要将其拟合成高斯曲线,然后从中获得半峰全宽。我可以做的 FWHM 部分,我已经有一个代码,但是我在编写适合高斯的代码时遇到了麻烦。

有没有人知道可以为我做这件事或能够为我指明正确方向的任何功能? (我可以对直线和多项式进行最小二乘拟合,但我不能让它适用于高斯)

如果它同时兼容 Octave 和 Matlab 也会很有帮助,因为我目前有 Octave,但要到下周才能访问 Matlab。

任何帮助将不胜感激!

【问题讨论】:

  • 你有一个峰值(只有 1 个高斯)吗?还是多个峰(多个重叠的高斯)?
  • 每个文件只有一个峰值。
  • 如果它只是一个峰值,取数字的平均值和标准差,这定义了您的样本正态分布。你试过吗?否则,如果您有统计工具箱,请使用 normfit()。
  • @Justin:你的第一句话是错误的。如果我在x=-100 周围有一堆数据点,y-values 对应于那里的标准法线,并且在x=-2 周围有一堆相似的值,那么所有这些点的平均值显然 not 为零,标准差显然不会为一。

标签: matlab octave curve-fitting gaussian


【解决方案1】:

也许这就是您要找的东西?不确定兼容性: http://www.mathworks.com/matlabcentral/fileexchange/11733-gaussian-curve-fit

来自其文档:

[sigma,mu,A]=mygaussfit(x,y) 
[sigma,mu,A]=mygaussfit(x,y,h)

this function is doing fit to the function 
y=A * exp( -(x-mu)^2 / (2*sigma^2) )

the fitting is been done by a polyfit 
the lan of the data.

h is the threshold which is the fraction 
from the maximum y height that the data 
is been taken from. 
h should be a number between 0-1. 
if h have not been taken it is set to be 0.2 
as default.

【讨论】:

  • 更适合作为评论,不是吗?
  • 虽然此链接可能会回答问题,但最好在此处包含答案的基本部分并提供链接以供参考。如果链接页面发生更改,仅链接答案可能会失效。 - From Review
  • @ScottHoltzman 感谢您的提醒,我已包含相关描述。
【解决方案2】:

我发现 MATLAB 的“fit”函数很慢,并使用带有内联高斯函数的“lsqcurvefit”。这是为了拟合高斯函数,如果你只想将数据拟合到正态分布,请使用“normfit”。

检查一下

% % Generate synthetic data (for example) % % %

    nPoints = 200;  binSize = 1/nPoints ; 
    fauxMean = 47 ;fauxStd = 8;
    faux = fauxStd.*randn(1,nPoints) + fauxMean; % REPLACE WITH YOUR ACTUAL DATA
    xaxis = 1:length(faux) ;fauxData = histc(faux,xaxis);

    yourData = fauxData; % replace with your actual distribution
    xAxis = 1:length(yourData) ; 

    gausFun = @(hms,x) hms(1) .* exp (-(x-hms(2)).^2 ./ (2*hms(3)^2)) ; % Gaussian FUNCTION

% % Provide estimates for initial conditions (for lsqcurvefit) % % 

    height_est = max(fauxData)*rand ; mean_est = fauxMean*rand; std_est=fauxStd*rand;
    x0 = [height_est;mean_est; std_est]; % parameters need to be in a single variable

    options=optimset('Display','off'); % avoid pesky messages from lsqcurvefit (optional)
    [params]=lsqcurvefit(gausFun,x0,xAxis,yourData,[],[],options); % meat and potatoes

    lsq_mean = params(2); lsq_std = params(3) ; % what you want

% % % Plot data with fit % % % 
    myFit = gausFun(params,xAxis);
    figure;hold on;plot(xAxis,yourData./sum(yourData),'k');
    plot(xAxis,myFit./sum(myFit),'r','linewidth',3) % normalization optional
    xlabel('Value');ylabel('Probability');legend('Data','Fit')

【讨论】:

  • 警告:由于代码使用lsqcurvefit,因此需要优化工具箱。
【解决方案3】:

我有类似的问题。 这是谷歌上的第一个结果,这里链接的一些脚本让我的 matlab 崩溃了。

最后我发现herematlab 内置了拟合函数,也可以拟合高斯。

看起来像这样:

>> v=-30:30;
>> fit(v', exp(-v.^2)', 'gauss1')

ans = 

   General model Gauss1:
   ans(x) =  a1*exp(-((x-b1)/c1)^2)
   Coefficients (with 95% confidence bounds):
      a1 =           1  (1, 1)
      b1 =  -8.489e-17  (-3.638e-12, 3.638e-12)
      c1 =           1  (1, 1)

【讨论】:

  • 注意fit不是内置的;它是曲线拟合工具箱的一部分
【解决方案4】:

直接拟合单个一维高斯是一个非线性拟合问题。你会发现现成的实现here,或here,或here for 2D,或here(如果你有统计工具箱)(你听说过谷歌吗?:)

无论如何,可能有一个更简单的解决方案。如果您确定您的数据y 将被高斯很好地描述,并且在整个x 范围内合理分布,您可以线性化问题(这些是方程式,而不是陈述):

   y = 1/(σ·√(2π)) · exp( -½ ( (x-μ)/σ )² )
ln y = ln( 1/(σ·√(2π)) ) - ½ ( (x-μ)/σ )²
     = Px² + Qx + R         

换人的地方

P = -1/(2σ²)
Q = +2μ/(2σ²)    
R = ln( 1/(σ·√(2π)) ) - ½(μ/σ)²

已经制作好了。现在,用(这些是 Matlab 语句)求解线性系统 Ax=b

% design matrix for least squares fit
xdata = xdata(:);
A = [xdata.^2,  xdata,  ones(size(xdata))]; 

% log of your data 
b = log(y(:));                  

% least-squares solution for x
x = A\b;                    

你用这种方式找到的向量x等于

x == [P Q R]

然后您必须对其进行逆向工程以找到均值 μ 和标准差 σ:

mu    = -x(2)/x(1)/2;
sigma = sqrt( -1/2/x(1) );

您可以与x(3) == R 进行交叉检查(应该只有小的差异)。

【讨论】:

  • 非常感谢。我只能通过 google 找到第一个链接,这不适用于我的数据,但第二个链接很有效。也感谢您的解释/方程式。 :D
  • @user1806676:我没有尝试过线性化方法,但至少数学是正确的。你应该在那里做一些试验和验证。
  • 尝试了线性化方法。效果很好。
  • +1。用于线性化并在找到系数之后取反对数。良好的最小二乘误差解决方案!
猜你喜欢
  • 1970-01-01
  • 2015-04-05
  • 1970-01-01
  • 1970-01-01
  • 2019-01-17
  • 2015-10-10
  • 1970-01-01
  • 2013-11-19
  • 1970-01-01
相关资源
最近更新 更多