How can I solve this function below?

parameters：a b c

KNOWN:(x1,y1),(x2,y2)... ...(xn,yn) How can I figure out all the parameters by using MATLAB?

Answer 1

[the cyclist]

Your code was fine, but your variables were pretty badly scaled, so I think the fit was behaving badly numerically. Here is some code where I rescaled your time variable for the fit (and then rescaled it back to display):

TIME_SCALE = 1000;

time = [1992   1993   1994   1995   1996   1997   1998   1999   2000   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010   2011   2012 ];
time = time/TIME_SCALE;

Y6711 = [0.001549451  0.001406211  0.001811241  0.001661079  0.002009743  0.001669221  0.002511118  0.001794201  0.003740947  0.00441657  0.004941945  0.004628348  0.004379816  0.005177686  0.004069535  0.004633793  0.007670759  0.009623098  0.009783128  0.010096395  0.012705525];

beta=nlinfit(time,Y6711,@curvefun,[0.5 0.5 0.5]);

a = beta(1);
b = beta(2);
c = beta(3);

%test the model
xx = min(time):(1/TIME_SCALE):max(time);
yy = a./(1+exp(b.*time+c));

plot(TIME_SCALE*time,Y6711,'o',TIME_SCALE*xx,yy,'r')

You are right that the initial guess is subjective. It is best to put in something that is approximately the right magnitude and sign if you can. For example, you know your parameter b has to be negative. But MATLAB figured it out here.

[Hui]

Thank you so much! There is no doubt that you have give me great help! But I still have some problems.Cause I have so many data to do the nonlinear fitting,and I'm not sure the ranges of the parametters so I can not set a good initial guess.

Is there any good idea? Thank you again!

[the cyclist]

I found your problem to be a bit tricky, but that is probably because I am not really an expert (and it is not something I do a lot of). I did spend some time thinking about it, though, because I found it quite interesting. Here's what I found.

Your parameter a doesn't really help the fit much, and seems to cause lots of problems because it is barely independent from c. I can get a very good fit with the function

function f=curvefun(a,x)
  f = 1./(1+exp(a(1).*x+a(2)));
end

instead.

When I use that, then I could do the following to make a guess at the initial parameters. Take the reciprocal of each side:

(1/y) = 1 + exp(bx+c)

Given your range of y, the 1 is negligible at first order, so ignore it:

(1/y) = exp(bx+c) [approx]

Take the log:

log(1/y) ~= bx + c [approx]

You can use a linear regression to estimate b and c:

b_linear = regress(log(1./Y6711'),[ones(size(time')),time']);

If you look at the relationship between this equation and the original one, you can see that a good guess at the initial parameters would be

b0 = [b_linear(2) b_linear(1)];

Putting this all together, I used the following code:

TIME_SCALE = 1000;

time = [1992   1993   1994   1995   1996   1997   1998   1999   2000   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010   2011   2012 ];
time = time/TIME_SCALE;

Y6711 = [0.001549451  0.001406211  0.001811241  0.001661079  0.002009743  0.001669221  0.002511118  0.001794201  0.003740947  0.00441657  0.004941945  0.004628348  0.004379816  0.005177686  0.004069535  0.004633793  0.007670759  0.009623098  0.009783128  0.010096395  0.012705525];

b_linear = regress(log(1./Y6711'),[ones(size(time')),time']);
b0 = [b_linear(2) b_linear(1)];

beta=nlinfit(time,Y6711,@curvefun,b0);

b = beta(1);
c = beta(2);

%test the model
xx = min(time):(1/TIME_SCALE):max(time);
yy = 1./(1+exp(b.*time+c));

figure
plot(TIME_SCALE*time,Y6711,'o',TIME_SCALE*xx,yy,'r')

function f=curvefun(a,x)
  f = 1./(1+exp(a(1).*x+a(2)));
end

This gave me a fit that looks just as good as the original (with one fewer parameter).

[Hui]

Thank you for your help!

Your later suggestion(two parameters model) is what I have done at the beginning. Cause the two parameters model cannot simulate some of my situations, so I need to improve this model by adding a parameter(three parameters model).As you said the questions maybe a bit tricky.

Thank you so much and you did give me great help undoubtedly.What's more, it makes me feel I am not working alone!

Answer 2

[the cyclist]

If you have the Statistics Toolbox, you could use the nlinfit() function to solve this.

[Hui]

Thank you for your answer I have try your suggestion, and this is my code below:

--------------------------------------------------------------------code
  function f=curvefun(a,x)
  f = a(1)./(1+exp(a(2).*x+a(3)));
  end
--------------------------------------------------------------------
time=[1992   1993   1994   1995   1996   1997   1998   1999   2000   2001   2002   2003   2004   2005   2006   2007   2008   2009   2010   2011   2012 ];
Y6711=[0.001549451  0.001406211  0.001811241  0.001661079  0.002009743  0.001669221  0.002511118  0.001794201  0.003740947  0.00441657  0.004941945  0.004628348  0.004379816  0.005177686  0.004069535  0.004633793  0.007670759  0.009623098  0.009783128  0.010096395  0.012705525];

beta=nlinfit(time,Y6711,@curvefun,[0.5 0.5 0.5 ]);
a=beta(1),b=beta(2),c=beta(3)
%test the model
xx=min(time):max(time);
yy=a./(1+exp(b.*time+c));
plot(time,Y6711,'o',xx,yy,'r')
--------------------------------------------------------------------code

RESULT: a =0.0051, b =-6.6651e+07, c =3.4666e+06

But I don't thing this result is right? And one of the parameters in nlinfit is subjectively defined. example:[0.5 0.5 0.5] Could you tell me where is the problems? Thank you so much!