Fitting to a Gaussian using nlinfit

[Ben]

I have the number concentration per particle size for a range of sizes and for 24 instrument channels. The plot seems to follow a normal distribution and I would like to fit the data to one. My data is not in the form of a probability density function since the y-axis is not probability, but rather concentration. Thus the equation I would like to fit to is: y = H.*exp(-.5.*((dp-mu)/sigma)^2), where H is the height of the distribution.

Following the methods outlined from the "Curve Fitting and Distribution Fitting" tutorial (http://www.mathworks.com/products/statistics/examples.html?file=/products/demos/shipping/stats/cfitdfitdemo.html), I attempted to use nlinfit to determine the optimal parameters H, mu, sigma. The result is not a good fit, however. I'm getting some errors when I run the code, namely, "Warning: Rank deficient, rank = 2, tol = 5.617004e-06." This, along with the terrible fit, makes me think I'm using this incorrectly.

Does anyone know what's going on here?

Here is my code:
%data
N_data = [
15284177237
23286616170
26033896374
39848389636
95606590189
1.95004E+11
4.1652E+11
6.90045E+11
9.16205E+11
1.13311E+12
1.34397E+12
1.52057E+12
1.70126E+12
1.87655E+12
1.84008E+12
1.75936E+12
1.527E+12
1.09247E+12
6.12101E+11
2.26027E+11
49697188712
17322725331
10247909480
5456170587];

dp_data = [
1.05E-08
1.14E-08
1.22E-08
1.31E-08
1.41E-08
1.51E-08
1.62E-08
1.75E-08
1.88E-08
2.02E-08
2.18E-08
2.34E-08
2.52E-08
2.72E-08
2.93E-08
3.15E-08
3.40E-08
3.67E-08
3.96E-08
4.27E-08
4.61E-08
4.98E-08
5.38E-08
5.81E-08];

%interpolate between points to get initial guess for height of dist.
dp_space = logspace(log10(dp_star_meas(1)*0.75),log10(dp_star_meas(length(dp_star_meas))*1.25),1000);
N_data_intr = interp1(dp_data,N_data,dp_space,'spline',0);
initial_guesses(1) = peak(N_data_intr);

%calculate mean and std for initial guess values
for i = 1:length(N_data)
mu = mu + dp_space(i)*N_data_intr(i);
N_tot = N_tot + N_data_intr(i);
end
mu = mu/N_tot;
for i = n:m
f = f + N_data_intr(i)*(dp_space(i) - mu)^2;
end
std = sqrt(f/(N_tot));
init_guesses(2) = mu;
init_guesses(3) = std;

%define function for nlinfit
norm_func = @(p,dp) p(1).*exp(-.5.*((dp-p(2))/p(3)).^2);

[bestfit,resid]=nlinfit(dp_data,log(N_data),@(p,x) log(norm_func(p,x)),init_guesses);

here is an image of the final fit versus the data: http://s2.postimage.org/xeu65jqyx/nlinfit_graph.png

I appreciate the help

[Torsten]

"Ben" wrote in message <kckmdj$8uo$1...@newscl01ah.mathworks.com>...

You fit the log of your data against the log of your fit function.
This introduces a distortion to the estimates of your parameters.
Why don't you try the direct way
[bestfit,resid]=nlinfit(dp_data,N_data,norm_func,init_guesses);
?

Best wishes
Torsten.

[Ben]

"Torsten" wrote in message <kcmih3$o3$1...@newscl01ah.mathworks.com>...

> You fit the log of your data against the log of your fit function.
> This introduces a distortion to the estimates of your parameters.
> Why don't you try the direct way
> [bestfit,resid]=nlinfit(dp_data,N_data,norm_func,init_guesses);
> ?
>
> Best wishes
> Torsten.

The reason why I used the log of my function is because of the reasons stated in the tutorial I linked to above (mainly that my function is not allowed negative values). I tried your suggestion, however, and the fit does look better but I'm still getting the same error message. Is there a way to restrict the allowed values for the parameters of nonlinfit? If anyone has any better suggestions for fitting a gaussian curve to this type of data, please let me know. I've been struggling with this for a while now.

[Torsten]

"Ben" wrote in message <kcmjis$4us$1...@newscl01ah.mathworks.com>...

Scale your N_data and your fit function by a common factor (e.g. 1e10).
Then the error message should disappear.

I obtain quite a good fit for
p(1)=1.95311e12
p(2)=2.85163e-8
p(3)=7.16019e-9

Best wishes
Torsten.

[Tom Lane]

Ben, nlinfit is doing a couple of things that are difficult in your problem.
First, it is trying to compute derivatives numerically by perturbing the
parameters a bit. Second, it is trying to measure convergence using the
change in parameters from one iteration to another. This is hard because you
have one parameter of the order of 1e12, and others of the order of 1e-8. A
simple way to avoid the problem is to build these scale factors into your
model, and make sure the parameters being estimated are of comparable
magnitude. The following shows how to do this on either the original or log
scale.

subplot(2,1,1)
init_guesses = [1 3 1];
norm_func = @(p,dp) 2e12*p(1).*exp(-.5.*((dp-1e-8*p(2))/(p(3)*1e-8)).^2);
[bestfit,resid]=nlinfit(dp_data,log(N_data),@(p,x)
log(norm_func(p,x)),init_guesses);
xx = linspace(min(dp_data),max(dp_data));
plot(dp_data, log(N_data),'bo', xx,log(norm_func(bestfit,xx)),'r-',
xx,log(norm_func(init_guesses,xx)),'g-')

subplot(2,1,2)
[bestfit,resid]=nlinfit(dp_data,N_data,@(p,x) norm_func(p,x),init_guesses);
plot(dp_data, N_data,'bo', xx,norm_func(bestfit,xx),'r-',
xx,norm_func(init_guesses,xx),'g-')

Further hint: when you post code, make sure it runs. You had multiple
undefined variables in your code. That's why I supplied my own init_guesses,
bypassing that part of your code.

-- Tom

[Ben]

Thank you to both Tom and Torsten, it seems like the function is doing what I want it to do now. I think I'm going to stick with the "non-log" version of the fit as it seems to give better results. I'll just have to hope it never goes crazy and starts giving negative values since I don't have the optimization toolbox.

Also, I apologize about the lack of defined variables, Tom. This is my first time posting code to a forum (which I know isn't a valid excuse) and I will make sure my code runs before posting next time.

[Nathan Orloff]

This will probably do it.

function [fitresult, zfit, fiterr, zerr, resnorm] = fmgaussfit(xx,yy,zz)
%fmgaussfit(xx,yy,zz) Custom Guassian Fit algorithm using lsqcurvefit
% A 3D gaussian fitting routine with error propagation and uncertainties.

%% Condition the data
[xData, yData, zData] = prepareSurfaceData( xx, yy, zz );
xyData = {xData,yData};

%% Set up the startpoint
[amp, ind] = max(zData); % amp is the amplitude.
xo = xData(ind); % guess that it is at the maximum
yo = yData(ind); % guess that it is at the maximum
ang = 45; % angle in degrees.
sy = 1;
sx = 1;
zo = median(zData(:))-std(zData(:));
xmax = max(xData)+2;
ymax = max(yData)+2;
xmin = min(xData)-2;
ymin = min(yData)-2;

%% Set up fittype and options.
Lower = [0, 0, 0, 0, xmin, ymin, 0];
Upper = [Inf, 90, Inf, Inf, xmax, ymax, Inf]; % angles greater than 90 are redundant
StartPoint = [amp, ang, sx, sy, xo, yo, zo];%[amp, sx, sxy, sy, xo, yo, zo];

tols = 1e-16;
options = optimset('Algorithm','levenberg-marquardt',...
'Display','off',...
'MaxFunEvals',5e2,...
'MaxIter',5e2,...
'TolX',tols,...
'TolFun',tols,...
'TolCon',tols ,...
'UseParallel','always');

%% perform the fitting
[fitresult,resnorm,residual] = ...
lsqcurvefit(@gaussian2D,StartPoint,xyData,zData,Lower,Upper,options);
[fiterr, zfit, zerr] = gaussian2Duncert(fitresult,residual,jacobian,xyData);

end

function z = gaussian2D(par,xy)
z = par(7) + ...
par(1)*exp(-(((xy{1}-par(5)).*cosd(par(2))+(xy{2}-par(6)).*sind(par(2)))./par(3)).^2-...
((-(xy{1}-par(5)).*sind(par(2))+(xy{2}-par(6)).*cosd(par(2)))./par(4)).^2);
end

function [dpar,zf,dzf] = gaussian2Duncert(par,resid,xy)
J = guassian2DJacobian(par,xy);
parci = nlparci(par,resid,'Jacobian',J);
dpar = (diff(parci,[],2)./2)';
[zf,dzf] = nlpredci(@gaussian2D,xy,par,resid,'Jacobian',J);
end

function J = guassian2DJacobian(par,xy)
x = xy{1}; y = xy{2};
J(:,1) = exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2);
J(:,2) = -par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(3).^2 - (2.*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2);
J(:,3) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2)./par(3)^3;
J(:,4) = (2.*par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2)./par(4)^3;
J(:,5) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2 - (2.*sind(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2);
J(:,6) = par(1).*exp(- (cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))).^2./par(3).^2 - (cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))).^2./par(4).^2).*((2.*cosd(par(2)).*(cosd(par(2)).*(y - par(6)) - sind(par(2)).*(x - par(5))))./par(4).^2 + (2.*sind(par(2)).*(cosd(par(2)).*(x - par(5)) + sind(par(2)).*(y - par(6))))./par(3).^2);
J(:,7) = ones(size(x));
end