How can I find the values for Settting Start Points when fitting point to a curve?

Hello:

I'm trying to fit some points

NUMBERof = transpose([8 12 36 37 15 152 21 182 107 169 235 324 394 ...
    462 545 707 655 769 832 838 812 849 864 950 932 879 604 637 743 ...
    757 683 605 510 619 517 567 523 551 585 565 410 399 430 435 440 ...
    367 378]);
DAY = transpose(1:length(NUMBERof));

with a Gaussian model:

fittype('a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2)',...
    'dependent',{'y'},...
    'independent',{'x'},...
    'coefficients', {'a1','b1','c1','a2','b2','c2','a3','b3','c3'});
myprevof = fit(DAY, NUMBERof,...
    'a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2)',...
    'start',[175 20 1 550 15 5 550 30 15]);

As you can see I've choosed a set of starting points (and I've tried several more) for the fitting, but the curve I always get is very bad.

I'd like to know if there's a method to calculate the better (or at least a good one) value for Star Points, and what is "Start Point" meaning to understand better the fitting.

Thank you in advance for your help.

Sergio.

Answer 1

[John D'Errico]

First, you don't understand how to use fittype. fittype does not set something that fit can use, UNLESS YOU PASS IN THE RESULTS OF FITTYPE INTO FIT.

Second, your data really does not merit three terms.

Finally, you can just use the gauss2 model form anyway.

F = fittype('gauss2')
F = 
     General model Gauss2:
     F(a1,b1,c1,a2,b2,c2,x) =  a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)

Now, call fit. But what parameter start estimates make sense there? PLOT YOUR DATA. LOOK AT THE PLOT. THINK ABOUT WHAT YOU SEE. ALWAYS DO THESE THREE THINGS.

plot(DAY,NUMBERof,'o')
grid on
xlabel DAY
ylabel casualties

Now you clearly have a peak around day 23. How wide is it? The width of that main peak is probably around 25-30 days, if we guess where the 10% points on the peak are. Divide that width by 4, to give a decent scale parameter in this form. How high is the main peak? I'll guess 800 at the center.

The second term is off to the right. Say it peaks around day 40. The width is wider, maybe I'll guess a width of 60 days. (Again, divide by 4 to be the parameter estimate.) How high? Say 300.

We don't need perfect numbers here. How do we use fit now? PASS IN F!

mdl = fit(DAY,NUMBERof,F,'start',[ 800, 23, 7.5, 300, 40, 15])

mdl = 
     General model Gauss2:
     mdl(x) =  a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)
     Coefficients (with 95% confidence bounds):
       a1 =       641.3  (247.6, 1035)
       b1 =       20.62  (19.74, 21.5)
       c1 =       8.827  (6.169, 11.48)
       a2 =       522.1  (422.6, 621.7)
       b2 =       36.25  (28.74, 43.76)
       c2 =       17.25  (8.059, 26.44)

My initial estimates were not dead on, but reasonably close. The secondary peak is disappointingly higher than I would like to see. But then people are too often foolish.

plot(mdl)
hold on
plot(DAY,NUMBERof,'bo')
xlabel DAY
ylabel casualties
grid on

Again, the fit seems entirely reasonable. The noise around the curve seems pretty random at this point, with little apparent lack of fit. That tells me any attempt to estimate a THIRD gaussian peak from that data would be a complete joke. You would be kidding yourself if you tried that, without considerably more and considerably better data.

In order to choose intelligent estimates for starting points, you need to understand what each parameter in the model does, and how it impacts the shape of the curve. Sorry, but there is no magical way to do that without understanding the model you are using.

[Sergio Sanz]

Thank you for your detailed analysis.

I agree with you when you say that I have to understand the model and what each parameter does. That's why I posted the question.

If you know where I can learn that in detail I'd appreciate very much your help.

Thank you again.

[Sergio Sanz]

Hello again:

Your method works fine, but I don't understand all the settings yet, and how to set the parameters in other cases.

So I'd appreciate if you give me some advice how or where I can learn it (websites or documents).

Here's my graphic using your idea of using gauss2 fitting:

Thank you again in advance.

Sergio

[John D'Errico]

As I said, you set the parameters by understanding the fundamental shape of the terms in your model. You have two terms here. Your model is the sum of TWO gaussian terms. What do each of those parameters do in the corresponding term? You should not be using that model if you have not a clue what it means. Yes, I know. Your teacher probably told you to use it.

So what does a single Gaussian term look like?

a1*exp(-((x-b1)/c1)^2)

PLOT IT. Set the parameters to simple versions. I might suggest

a1 = 1;
b1 = 0;
c1 = 1;

Is this simple function symmetric around some point? Should it be?

With those parameters set as I suggest, where does that function attain its maximum? What happens when x == b1?

How does a change of b1 to some value other than zero impact the shape of the curve?

What is the value of that function when x==b1?

If you change the value of a1 from 1, how does that impact the maximum?

How does that function drop off as x --> inf? At what value does that shape drop off to say 10% of the maximum value? How is that a function of c1? That is, if c1 is larger or smaller than 1, how would the function change its shape? Does the shape really change, or is c1 just a scaling parameter?

The point is, ALL of those parameters do not really change the fundamental shape of that single gaussian term. They just scale it in some way, or translate it along the x axis.

Next think about ANY function. Suppose you have some other function. I'll just call it F(x). I really don't care what the shape is at this point.

What does a1*F(x) look like? How does that compare to F(x)?

Next, what does F(x - b1) look like? How does that compare to F(x)?

Next, what does F(x/c1) look like? How does THAT compare to F(x)?

Next, put it all together. If you understand what I asked in the last three questions, then what does

a1*F((x - b1)/c1)

look like, in comparison to F(x)?

If you understand what I have said above, you don't need to read websites or documents. If not, then go back and read what I have said again. Think about what each of those parameters mean, about how they impact the shape of the function they are applied to.

[Sergio Sanz]

Hello:

Great way to understand the method!

You're a great teacher!

Thank you very much!

:-)

[Sergio Sanz]

Hello again:

It's easy to understand the behaviour of the curve when changing the coefficients, but my troubles are about the STARTING VALUES to plot the graphics in Matlab.

I can't see the relationship between the two things.

Sorry for questioning again.

[John D'Errico]

Here I'm a bit confused as to what you are asking. You only need to supply starting values to perform the fit. Essentially, the fitting routine needs to know a guess as to reasonable set of parmeters for the model. Then it refines those estimates. But you don't need starting values for the plot.

For example, here is the model I fit to your data.

F = fittype('gauss2');
mdl = fit(DAY,NUMBERof,F,'start',[ 800, 23, 7.5, 300, 40, 15])
mdl = 
     General model Gauss2:
     mdl(x) =  a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)
     Coefficients (with 95% confidence bounds):
       a1 =       641.3  (247.6, 1035)
       b1 =       20.62  (19.74, 21.5)
       c1 =       8.827  (6.169, 11.48)
       a2 =       522.1  (422.6, 621.7)
       b2 =       36.25  (28.74, 43.76)
       c2 =       17.25  (8.059, 26.44)

plot(mdl)
hold on
plot(DAY,NUMBERof,'bo')
xlabel DAY
ylabel casualties
grid on

As you can see, I supplied starting values. But there were needed for the fit only. After the fit is performed, that model structure contains the fitted parameters. Here we can see it does contain the fitted values:

mdl.a1
ans =
       641.28

At this point, no starting values are needed. You can do the plot directly from the model.

So are you talking about the prediction intervals around the model as you plottted it? The predint code does allow you to supply some additional information to it, but the default values will be reasonabe to use as they are. I could explain some of them, if that is your question. I'm just not sure what you mean about starting values for the plot.