How to find the initial parameters for nlinfit?

I have such set data:

t =
   [ 0
  0.0058
  0.0116
  0.0174
  0.0233
  0.0291
  0.0349
  0.0407
  0.0465
  0.0523
  0.0581
  0.0640
  0.0698
  0.0756
  0.0814
  0.0872]
y=1.0e-04. *
  [0.0177
  0.0977
  0.1102
  0.0797
  0.0174
  0.0010
  0.0003
  0.0003
  0.0001
  0.0000
  0.0000
  0.0000
  0.0000
  0.0000
  0.0000
  0.0000]

I want to use scaled gamma distribution function to fit them. The problem is I do not know how to set the suitable initial guess for parameters. I tried many sets but all failed when calling the nlinfit function. The code is just like this:

lambda0=[0.5 1 1]  %this set is not right

lambda2 = nlinfit(t,y,@low_pass_fun,lambda0);

% % show the filter
figure(22);
SetGrayscale;
plot(t,y,'k*');hold on;
plot(t,low_pass_fun(lambda2,t),'g');hold off;

function yfun = low_pass_fun(lambda,t)
    esp=1e-30;
    yfun = lambda(1)*(t+esp).^(lambda(2)-1) .* exp(-t./lambda(3))./(gamma(lambda(2))*lambda(3)^lambda(2));  % scaled gamma distribution

end

How to find the initial parameters? Thanks a lot.

Answer 1

[Tom Lane]

I'm writing this as a separate answer just so I can use the mark-up feature. I took the code you supplied, including the data to exactly the precision you showed in your question, and I ran this:

lambda0 = [2e-7,1.3,.009]
lambda2 = nlinfit(t,y,@low_pass_fun,lambda0)

I got no warnings. I'm surprised you did. Maybe it could be related to the OS -- I'm not sure.

Warnings early in the fitting process are less of a concern. The nlinfit function may get itself out of a bad situation as the iterations proceed. You could try taking the lambda2 values, then run nlinfit again using those as the starting values. If you get no warnings from that, you might have more confidence that things are okay.

Answer 2

[Tom Lane]

Here's how I approach problems like this. Other variants may work for you. I'll approximate the mean as the sample mode (location of peak), and the standard deviation is something like the width of the peak or 0.01. I can look up the gamma mean and standard deviation as a function of the two parameters a and b. Then I can solve for a and b.

>> t(find(y==max(y)))              % mode
ans =
    0.0116
>> % a*b = .0116, sqrt(a)*b = .01  % known formula for mean and std dev
>> a = (.0116/.01)^2               % solve for a
a =
    1.3456
>> b = .0116/a                     % solve for b
b =
    0.0086

Now I need the scale factor, which is the integral, or roughly:

>> sum(y(2:end).*diff(t))
ans =
   1.7806e-07

So my starting values could be something like

lambda0 = [2e-7,1.3,.009]

[Yunhui]

Thank you so much.

I tried your method and received a lot of warnings of nlinfit:

Warning: Rank deficient, rank = 2, tol = 3.0645e-04.

> In nlinfit>LMfit at 294

In nlinfit at 166

..........................

But the figure shows the fitting is quite good.

Why this time rank deficient does not affect the fitting? But when I met such kind of warnings, the fitting always fails.

Thank you again.

[Tom Lane]

I would have to see exactly what you did. I tried your problem, including your data, using the lambda0 that I calculated. No warnings. It is possible to get a good fit with warnings; for example if you had redundant parameters that could happen.

[Yunhui]

I am using Matlab R2010a. I tried many times, using that data and the lambda0 that you calculated, but always received many warnings.

I have downloaded R2012a also. But when I run this code, it said : Undefined function 'nlinfit' for input arguments of type 'function_handle'.

And I found nlinfit.m is in different folder compared to R2010a.

In Matlab R2010a, nlinfit.m is in '/Applications/MATLAB_R2010a.app/toolbox/stats' (I am using mac system)

while in Matlab R2012a, nlinfit.m is in '/Applications/MATLAB_R2012a.app/toolbox/stats/stats' .

How can I make it run in R2012a?

Thanks a lot.

[JAY R]

nice description sir.

Answer 3

[Tom Lane]

I just tried your problem in R2010a using your code and my lambda0, and all worked well. In R2012a the directory structure of the Statistics Toolbox is indeed different. You might want to try "which pathdef" and delete the file you find if it sets the path incorrectly. You could also "cd(matlabroot)" and "cd toolbox/stats/stats" if you just want to try it once, but I do not recommend this as a long-term solution.

[Yunhui]

I use your method to fix the pathdef file. Now I can run R2012a version! Thanks a lot.

But still receive a lot of warnings. What is my wrong?

Is it because of the different operating system? But anyway, thank you so much for your detailed and patient answers!

[Yunhui]

Just now I asked my friend to try it. She is using Mac and Matlab R2011b. She got the same set of warnings. And I also tried on a windows machine( win 7) and Matlab R2012a, the same thing happened.

So strange that you got no warnings.

And I tried using lambda2 as starting value to ran nlinfit again, still a lot of warnings.

But the figure shows the fitting is good. Thanks a lot for you to teach me such a good way to find initial paramters.

[JAY R]

Good Morning sir, This is regarding your answer I have read http://in.mathworks.com/matlabcentral/answers/35074-how-to-find-the-initial-parameters-for-nlinfit. you have explain

>> t(find(y==max(y))) % mode ans = 0.0116

>> % a*b = .0116, sqrt(a)*b = .01 % known formula for mean and std dev

>> a = (.0116/.01)^2 % solve for a a = 1.3456 >> b = .0116/a % solve for b b = 0.0086 Now I need the scale factor, which is the integral, or roughly:

>> sum(y(2:end).*diff(t)) ans = 1.7806e-07 So my starting values could be something like

lambda0 = [2e-7,1.3,.009]

***********************************

I have a small doubt? If I have one variable let us say vec1 is 32x5 and vec2 is 32x1 then how to find the initial guesses if the function is of five paramter ?

You may take random value for vec1 and vec2 for the explanation. Thanks in advance

[Tom Lane]

The material you quote uses specific knowledge of the function to be fit, namely the form of the gamma distribution, to figure out starting values. The initial guesses in your case would depend on the form of the function you intend to fit.