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

Question

Doug Barrett 2013-11-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/107964-using-nlinfit-to-fit-a-gaussian-pdf-to-x-y-paired-data

评论： Doug Barrett 2013-12-2

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');

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Image Analyst 2013-11-28

编辑：Image Analyst 2013-11-28

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 2013-11-29

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/107964-using-nlinfit-to-fit-a-gaussian-pdf-to-x-y-paired-data#answer_116721

编辑：Matt J 2013-11-29

在 MATLAB Online 中打开

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,

3 个评论
显示 1更早的评论隐藏 1更早的评论

Doug Barrett 2013-11-29

编辑：Doug Barrett 2013-11-29

在 MATLAB Online 中打开

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 2013-11-29

编辑：Matt J 2013-11-29

在 MATLAB Online 中打开

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 2013-12-2

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.

请先登录，再进行评论。