Need help in using nlinft.

I'm having trouble in fitting a logistic function using 'nlinfit'.

Code:

x = [1912 1921 1930 1936 1956 1960 1968 1983 1988 1991 1994 1996 1999 2005 2007]';
y = [10.6000 10.4000 10.3000 10.2000 10.1000 10.0000 9.9500 9.9300 9.9200 9.8600 9.8500 9.8400 9.7900 9.7700 9.7400]';

fun = @(p,x)p(1)+p(2).*exp(-p(3).*x);     % logistic function

p0 = [1 1 1];      % initials

p = nlinfit(x, y, fun, p0);

xx = 1900:10:2200;
yy = feval(fun, p, xx);

plot(x,y,'o',xx,yy)

The result is horrible. I have tried testing many initials, but the fitting curves were not good at all.

Is it because the data set is too small? But those are the only data I can acquire.

Thank you so much for your help!

Answer 1

[the cyclist]

Your variables are badly scaled, so are going to be numerically unstable when you put them in the logit function. Try rescaling your x variable:

XSCALE = 1000;

x = [1912 1921 1930 1936 1956 1960 1968 1983 1988 1991 1994 1996 1999 2005 2007]';
x = x/XSCALE;

y = [10.6000 10.4000 10.3000 10.2000 10.1000 10.0000 9.9500 9.9300 9.9200 9.8600 9.8500 9.8400 9.7900 9.7700 9.7400]';

fun = @(p,x)p(1)+p(2).*exp(-p(3).*x);     % logistic function

p0 = [1 1 1];      % initials
p = nlinfit(x, y, fun, p0);

xx = 1900:10:2200;
xx = xx/XSCALE;
yy = feval(fun, p, xx);

figure
plot(XSCALE*x,y,'o',XSCALE*xx,yy)

[mnt]

Thx! But the fit is not satisfying...

[the cyclist]

Here is another refinement of the code:

XSCALE = 1000;
YSCALE = 10;

x = [1912 1921 1930 1936 1956 1960 1968 1983 1988 1991 1994 1996 1999 2005 2007]';
x = x/XSCALE;

y = [10.6000 10.4000 10.3000 10.2000 10.1000 10.0000 9.9500 9.9300 9.9200 9.8600 9.8500 9.8400 9.7900 9.7700 9.7400]';
y = y/YSCALE;

fun = @(p,x)p(1)+exp(p(2)-p(3).*x);     % logistic function

p0 = [1 30 15];      % initials
p = nlinfit(x, y, fun, p0);

xx = 1900:10:2020;
xx = xx/XSCALE;
yy = feval(fun, p, xx);

figure
plot(XSCALE*x,YSCALE*y,'o',XSCALE*xx,YSCALE*yy)

You'll notice I did a couple things here. First, I also scaled your y variable. This helped a little. But the code was still giving a warning that it exceeded the maximum number of iterations.

You can set the value of that maximum number in the code, but instead of doing that, I took the estimates of p from the fit, and set those values to be the new initial estimates (approximately). That's why you'll see that I have set

p0=[1 30 15]

in the code.

The fit is much better now.

[mnt]

Thank you so much! Learned a lot from you!

[the cyclist]

The best form of thank you is to accept the answer, which may also help future users.