Using nlinfit to fit a Gaussian pdf to x,y paired data

I am trying to use Matlab's nlinfit function to estimate the best fitting Gaussian pdf for x,y paired data. In this case, x is a range of 2D orientations and y is the probability of a "yes" response in a yes-no task.

I have borrowed the form for a Gaussian (@norm_funct) from relevant posts and I'd like to estimate the best fitting pdf to obtain its mean and standard deviation. At the moment, the fitted function appears to be incorrectly scaled and is not smooth (see plots)

The figure below, which was produced using Prism to fit 2 pdfs, is what I'd like to replicate in Matlab - many thanks

x = -30:5:30;
y = [0,0.20,0.05,0.15,0.65,0.85,0.88,0.80,0.55,0.20,0.05,0,0;];

% plot raw data
figure(1)
plot(x, y, ':rs');
axis([-35 35 0 1]);

% initial paramter guesses (based on plot)
initGuess(1) = max(y);  % amplitude
initGuess(2) = 0;       % mean centred on 0 degrees
initGuess(3) = 10;      % SD in degrees

% equation for Gaussian distribution
norm_func = @(p,x) p(1) .* exp(-((x - p(2))/p(3)).^2);

% use nlinfit to fit Gaussian using Least Squares
[bestfit,resid]=nlinfit(y, x, norm_func, initGuess);

% plot function
xFine = linspace(-30,30,100);
figure(2)
plot(x, y, 'ro', x, norm_func(xFine, y), '-b');

[Image Analyst]

Should I delete your apparent duplicate http://www.mathworks.com/matlabcentral/answers/107961-using-nlinfit-to-fit-a-gaussian-pdf-to-x-y-paired-data? Or are they really different questions? Why don't we have just one thread so you don't have to look in two different places for answers?

Answer 1

[Matt J]

Instead of

 plot(x, y, 'ro', x, norm_func(xFine, y), '-b');

I think you want

 plot(x, y, 'ro', xFine, norm_func(bestfit,xFine), '-b');

From this, I obtain the following (but from fminsearch, not nlinfit), and it looks fine to me,

[Doug Barrett]

Thanks Matt, looks just-the-job!

Could I ask how you implemented fminsearch for norm_func?

I have tried the code below but I get an error meassage indicating I'm not entering enough parameters. I'm calling 'norm_func' from another file rather than anonymously and I'd like to minimise the (sum of squares) difference between the predicted and observed data (f(x) and y).

% non linear fit for x,y data where x = degrees of rotation and y =  
% probability of a yes response

clear all

x = -30:5:30;
y = [0,0.20,0.05,0.15,0.65,0.85,0.88,0.80,0.55,0.20,0.05,0,0;];

% plot raw data
figure(1)
plot(x, y, ':rs');
hold on

% initial paramter guesses (based on plot)
initGuess(1) = max(y);  % amplitude
initGuess(2) = 0;       % mean centred on 0 degrees
initGuess(3) = 10;      % SD in degrees

options = optimset('fminsearch');
bestfit = fminsearch('norm_func', initGuess, options, [x,y]);

% plot function
xFine = linspace(-30,30,100);
plot(xFine, norm_func(bestfit,xFine), '-b');
axis([-35 35 0 1]);

Here's norm_func

function SS = norm_func(initGuess, x, y)

% Gaussian pdf
pPred = initGuess(1).* exp(-((x - initGuess(2))/initGuess(3)).^2);
pObsv = y;
SS = sum((pObsv - pPred).^2);  % LS diff

end

Thanks!

[Matt J]

Doug,

Here's how I modified your original code to work with fminsearch. However, I wouldn't conclude that fminsearch is preferable to nlinfit. Fminsearch is, in fact, a rather ad hoc and less robust solver. The only reason I used it instead of nlinfit is because I don't have the Stats Toolbox.

x = -30:5:30;
y = [0,0.20,0.05,0.15,0.65,0.85,0.88,0.80,0.55,0.20,0.05,0,0;];

% plot raw data
figure(1)
plot(x, y, ':rs');
axis([-35 35 0 1]);
% initial paramter guesses (based on plot)
initGuess(1) = .9*max(y);  % amplitude
initGuess(2) = 0;       % mean centred on 0 degrees
initGuess(3) = 10;      % SD in degrees
% equation for Gaussian distribution
norm_func = @(p,x) p(1) .* exp(-((x - p(2))/p(3)).^2);
% use nlinfit to fit Gaussian using Least Squares
[bestfit,resid]=fminsearch(@(p) norm(norm_func(p,x)-y), initGuess);
% plot function
xFine = linspace(-30,30,1000);
figure(2)
plot(x, y, 'ro', xFine, norm_func(bestfit,xFine), '-b');

[Doug Barrett]

Hi Matt,

Thanks for that, it's very neat. I'll look into the nlinfit more now I have a better idea of how to specify and call the function - I might look at 'mle' too (statstoolbox), as it looks as though I could use a binomial function to fit the accracy data and specify upper and lower bounds for the estimates.

All the best,

Doug.