Issues with fminsearch

Hello everyone,

To start my question let me just focus on the fact that I'm a novice when it comes to programming and especially Matlab.

I'm currently struggling in finding a way to non-linear fit my data using two Gaussian plots. I would like to use fminsearch to achieve this task. Right now this is the code I was able to generate. First I generate a gaussian curve with some noise then would like to model it.

function y = gauss(x)

for i = 1:101
g(i) = 100*exp ( - ((i - 50)^2) / (2*20^2))
end

y = g + randn(1,101);

subplot(1,2,1)
plot(g)
subplot(1,2,2)
plot(y)

objfun = @(x) 50*exp( -((x-70)^2) / (2*10^2) );

[x fval] = fminsearch(objfun, y);

end

Thanks again for all the help!

[SCd merging] Here is what I'm trying to achieve.

1) I created a Gaussian peak containing random noise. This could be generated using this code.

for i = 1:101 
g(i) = 100*exp ( - ((i - 50)^2) / (2*20^2)) 
end 
y = g + randn(1,101);

2) I would like to non-linear fit this Gaussian peak using fminsearch. I want to minimize the sum of square separating my data set from the model. The model I would like to implement is a Gaussian.

3) At page 33 of this document http://perso.unige.ch/~manuel/Matlab/MatLab1.pdf, there is a figure describing exactly what I want to do.

[/merge]

[SCd merging]

However, here is the bigger picture on what I'm trying to do.

I would like to fit 2 Gaussians on this 2d plot: http://imageshack.us/photo/my-images/594/2dscan.jpg/ using the program you just wrote. This data set is 1:101 matrix.

[/merge]

Answer 1

[Walter Roberson]

Change

objfun = @(x) 50*exp( -((x-70)^2) / (2*10^2) );

to

objfun = @(x) 50*exp( -((x-70).^2) / (2*10^2) );

[Sean de Wolski]

Philippe writes:

Tried it. Unfortunately I have the following error out.

Subscripted assignment dimension mismatch.

Error in fminsearch (line 191) fv(:,1) = funfcn(x,varargin{:});

Error in gauss (line 24) [x fval] = fminsearch(objfun, y)

Answer 2

[Sean de Wolski]

Edit

So you have a vector, 1:101. You do a calculation of this with a function, g, and add some noise to the result. What now do you want to calculate? You know everything!

The thing you would be trying to calculate would be parameters (coefficients) used in the function g (e.g. the 100 or the 50 if those has been unknowns and you had only been given x/y and the function.). But since there are no unknowns, this doesn't make sense.

For the purpose of explaining this to you, consider this change to your problem statement: given a vector vec of inputs, a function fun with some unknowns (where the 50 and the 100 were), and a set of noisy data g of that function evaluated at vec + noise. We now want to figure out those parameters:

vec = 1:101;
fun = @(x)100*exp ( - ((x - 50).^2) ./ (2*20^2)); %known evaluating
g = fun(vec);
gnoisy = g+randn(size(vec));

fun2 = @(a,x)a(1)*exp ( - ((x - a(2)).^2) ./ (2*20^2)); %unknowns
nlinfit(vec,gnoisy,fun2,[105 48]) 
%inputs vec, noisy output, function where we don't know a1, a2, initial guess

[Sean de Wolski]

Philippe wrote:

I got this message error when running your program:

Error using nlinfit (line 120) Error evaluating model function '@(x)100*exp(-((x-50).^2)./(2*20^2))'.

Caused by: Error using @(x)100*exp(-((x-50).^2)./(2*20^2)) Too many input arguments.

nlinfit is a Nonlinear least-squares regression of a data set. However, fminsearch uses the Nelder-Mead simplex method, which is what I would like to use.

[Sean de Wolski]

That was a copy error on my part, should've been fun2. It's fixed.

It's still not clear what you're trying to do and thus not at all apparent why yo want fminsearch.

To do a curve fit with this data, nlinfit is a good tool and is designed for this type of problem.

Also, please post replies in comments not answers.

[Philippe]

What is the [105 48] for? Is it possible to plot the final model and compare it with the initial data?

Based on your corrected program, I would like to model this 2d plot : http://imageshack.us/photo/my-images/594/2dscan.jpg/. This data set is 1:101 matrix which I would like to fit to two Gaussian curves.

[Sean de Wolski]

The 105 48 is an initial guess for the paremeters. As you can see the result is pretty close to the actual used [100 50]. Any deviations account for the random addition.

To overlay the plots, call fun2 with the original x (vec) and the output from nlinfit (a)

a = nlinfit(...

g2 = fun2(a,vec)

[Philippe]

OK! COOL! Thanks for that one!

However, here is the bigger picture on what I'm trying to do.

I would like to fit 2 Gaussians on this 2d plot: http://imageshack.us/photo/my-images/594/2dscan.jpg/ using the program you just wrote. This data set is 1:101 matrix.

Answer 3

[Philippe]

Based on your reply and my programming knowledge (which is very limited), I'm not sure how I can achieve this...

[Sean de Wolski]

Explain, in English with no code, what you're trying to do. What do you want to find? What is unknown? What is known? What are your inputs? Etc. Then supplement that with code or pseudocode if necessary.

If you can communicate the goals to us, we can help you with the syntax/implementation.

[Sean de Wolski]

Philippe writes:

If you superpose y with x on the same plot, you will see that it is not well fitted...

[Sean de Wolski]

You need to prep you objfun to be something to be minimized. Right now it will be a minimum of 0 at any large value (abs(x)>a few hundred). It is improbable if not impossible that it would converge to the desired solution.

objfun(inf),objfun(800),objfun(-723)

You might want to try nlinfit() (in the statistics toolbox). It will set this type of fit up for you.