曲線擬合——（1）高斯曲線

作者：桂

時間：2017-03-13  21:23:57

鏈接：http://www.cnblogs.com/xingshansi/p/6545162.html

---

前言

 本文主要是上一篇文章的補充，主要針對常用正態分布曲線擬合，文中內容多有借鑒他人，最后一並給出鏈接。

 

一、理論分析

對於正態分布：

$f(x) = \frac{1}{{\sqrt {2\pi } \sigma }}{e^{ - \frac{{{{(x - \mu )}^2}}}{{2{\sigma ^2}}}}}$

假設數據點{$x_i$,$y_i$}（$i = 1,2,3,...N$）符合正態分布曲線，對其進行擬合（曲線擬合不同於分布擬合，需要乘以幅度$A$）,給出准則函數：

對准則函數$J_0$求解即可實現參數估計。

由於求導比較復雜（可以借助Mathmatica/Maple），因此這里換一個思路：如果$e^x$—>$y$，則$x$—>$lny$，重新定義准則函數：

此時，變成了{$x_i$,$ln(y_i)$}的一元二次最小二乘擬合問題（此步簡便，直接借助MATLAB的polyfit，不再細說）。假設擬合結果為：

對應參數估計：

二、代碼實現

　　A-編程

根據上文理論分析，直接給出代碼：

clc;clear all;close all;
%generate the orignal data
mu = 0;
sig2 = 2;
A = 4;
x=-10:.1:10;
y=A/sqrt(2*pi*sig2)*exp(-(x-mu).^2/sig2/2)+.05*randn(1,length(x));
subplot 211
scatter(x,y,'k');grid on;

%%curve fitting
%1-Gauss distribution
p = polyfit(x,log(y),2);
sig2_est = -1/2/p(1);
mu_est = sig2*p(2);
A_est = exp(mu^2/2/sig2+p(3))*sqrt(2*pi*sig2);

%plot
subplot 212
scatter(x,y,'k');hold on;
grid on;
plot(x,A/sqrt(2*pi*sig2)*exp(-(x-mu_est).^2/2/sig2),'r--','linewidth',2);

結果圖：

　　B-自帶函數

如果單從實現角度，可以直接調用：

clc;clear all;close all;
%generate the orignal data
mu = 0;
sig2 = 2;
A = 4;
x=-10:.1:10;
y=A/sqrt(2*pi*sig2)*exp(-(x-mu).^2/sig2/2)+0.05*randn(1,length(x));
subplot 211
scatter(x,y,'k');grid on;

%%curve fitting
%1-Gauss distribution
f = fittype('a*exp(-((x-b)/c)^2)');
[cfun,gof] = fit(x(:),y(:),f);
yo = cfun.a*exp(-((x-cfun.b)/cfun.c).^2);
%plot
subplot 212
scatter(x,y,'k');hold on;
grid on;
plot(x,yo,'r--','linewidth',2);

對應結果圖：可見二者都可以實現擬合。

 

三、擬合優化

推導時，我們用了一個前提：如果$e^x$—>$y$，則$x$—>$lny$。對於接近於0處的噪聲，$lny$顯然會將噪聲放大。現在增加噪聲，觀察編寫的code與自帶兩種結果：

可以看到，對於$y_i$接近0處的噪聲，$lny$會將噪聲放大。擬合結果非常不理想，現在考慮對編程擬合方法進行優化：

既然在映射到$ln$函數時，$y_i$接近於0點處的噪聲會放大，考慮只通過數值較大的點進行擬合，而直接舍去接近0的點。即：添加閾值Th，僅考慮部分樣本點。

給出優化后的代碼：

clc;clear all;close all;
%generate the orignal data
mu = 3;
sig2 = 2;
A = 40;
xold=-10:.1:10;
yold=A/sqrt(2*pi*sig2)*exp(-(xold-mu).^2/sig2/2)+0.5*randn(1,length(xold));
subplot 211
scatter(xold,yold,'k');grid on;
%%curve fitting
%select data
x = [];y = [];
Th = 0.3;%Threshold
for n = 1:length(xold)
    if yold(n)>max(yold)*Th;
        x = [x,xold(n)];
        y = [y,yold(n)];
    end
end

%1-Gauss distribution
p = polyfit(x,log(y),2);
sig2_est = -1/2/p(1);
mu_est = sig2*p(2);
A_est = exp(mu^2/2/sig2+p(3))*sqrt(2*pi*sig2);

%plot
subplot 212
scatter(xold,yold,'k');hold on;
grid on;
plot(xold,A/sqrt(2*pi*sig2)*exp(-(xold-mu_est).^2/2/sig2),'r--','linewidth',2);

優化結果：　

問題得到改善。

參考：

高斯擬合：http://www.cnblogs.com/linkr/p/3632032.html