How can I solve this function below?

1 view (last 30 days)
Hui
Hui on 26 Apr 2014
Commented: Hui on 30 Apr 2014
parameters:a b c
KNOWN:(x1,y1),(x2,y2)... ...(xn,yn) How can I figure out all the parameters by using MATLAB?

Accepted Answer

the cyclist
the cyclist on 27 Apr 2014
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.
  3 Comments
the cyclist
the cyclist on 29 Apr 2014
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
Hui on 30 Apr 2014
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!

Sign in to comment.

More Answers (1)

the cyclist
the cyclist on 26 Apr 2014
If you have the Statistics Toolbox, you could use the nlinfit() function to solve this.
  1 Comment
Hui
Hui on 27 Apr 2014
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!

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!