Re: Dancing "a la Levenberg-Marquardt" to get the best Logistic Model.
- To: mathgroup at smc.vnet.net
- Subject: [mg124848] Re: Dancing "a la Levenberg-Marquardt" to get the best Logistic Model.
- From: Gilmar Rodriguez-pierluissi <peacenova at yahoo.com>
- Date: Thu, 9 Feb 2012 05:39:45 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201202081033.FAA06594@smc.vnet.net> <4F3298B0.5080909@wolfram.com>
- Reply-to: Gilmar Rodriguez-pierluissi <peacenova at yahoo.com>
Thank you for your valuable help Dr. Glosemeyer!
Sincerely,
Gilmar Rodriguez Pierluissi
From: Darren Glosemeyer <darreng at wolfram.com>
To: Gilmar Rodriguez-pierluissi <peacenova at yahoo.com>
Cc: mathgroup at smc.vnet.net
Sent: Wednesday, February 8, 2012 10:45 AM
Subject: [mg124848] Re: [mg124829] Dancing "a la Levenberg-Marquardt" to get the best Logistic Model.
On 2/8/2012 4:33 AM, Gilmar Rodriguez-pierluissi wrote:
> Dear Math Group:
> (**I start with the following population values between 1972 to 2008:**)
> In[2]:= popvalues = {217928, 219129, 221577, 227481, 231748, 233514,
> 232857, 233664, 235228, 240526, 243310, 249587, 250128, 253383,
> 257751, 261999, 258229, 262567, 263272, 267643, 272468, 274035,
> 276154, 278323, 282606, 289505, 295243, 293956, 294410, 296399,
> 297382, 298289, 299248, 296785, 299359, 300184, 299993};
> In[3]:= L = Length[popvalues];
> (** The highest (projected) value that the population can reach is: **)
> In[4]:= pop2020 = 304909;
> (** Assemble the data to do build a "scatter plot" **)
> In[5]:= popvalues2D =
> Join[Table[{i, popvalues[[i]]}, {i, 1, L}], {{49, pop2020}}];
> In[6]:= plt1 = ListPlot[popvalues2D]
> Out[6]= (** Plot ommited **)
> (** Let: **)
> In[7]:= K = pop2020
> Out[7]= 304909
> In[8]:=
> Subscript[P, 1] = popvalues[[1]]
> Out[8]= 217928
> (** I'm attempting to use a Population Logistic model similar to one \
> found in (where else?) Wikipedia:
> http://en.wikipedia.org/wiki/Logistic_function under the title: "In \
> ecology: modeling population growth". **)
> (** Since I need this model to satisfy Logistic[1]= Subscript[P, 1] \
> and Lim t -> Infinity Logistic[t] = K; I came up with the following \
> version of the Logistic Model to handle the above data set \
> appropriately: **)
> (** Logistic[t_]=(K Subscript[P, 1]Exp[rt])/(K Exp[r]+ Subscript[P, \
> 1](Exp[er]-Exp[r])); **)
> (** If you inspect this model ("by hand") you will see that \
> Logistic[1]= Subscript[P, 1] (the first population data point). Using \
> L'Hopital's Rule; one can show that Lim t -> Infinity (Logistic[t]) = \
> K; by taking the derivative of the numerator and denominator with \
> respect to t and performing the appropriate cancellations. Again; K \
> is the highest value that the population can reach "by design". **)
> (** Logistic[t_]=(Subscript[P, 1] E^(r*t))/(E^r+ Subscript[P, 1] \
> (E^(r*t)- E^r)/K); **)
> (** The model is equivalent to: **)\[AliasDelimiter]
> In[13]:= Logistic[t_] = ( Subscript[P, 1] Exp[r t])/(
> Exp[r] + Subscript[P, 1] (Exp[r t] - Exp[r])/K);
> (** I' m expecting Logistic[1] = 217928 and indeed : )
> In[14]:= Logistic[1]
> Out[14]= 217928
> (** but, unfortunately; **)
> In[16]:= Limit[Logistic[t], t -> Infinity]
> Out[16]= Limit[(217928 E^(r t))/(
> E^r + (217928 (-E^r + E^(r t)))/304909), t -> \[Infinity]]
> (** and: **)
> In[15]:= Logistic[49]
> Out[15]=
> = (217928 E^(49 r))/(E^r + (217928 (-E^r + E^(49 r)))/304909)
> (** I can see the the function Logistic[t] requires to be "herded" \
> (somehow) so that cancellations of terms can take place. Perhaps \
> using "Hold[]" and "ReleaseHold[]; I just don't know how. **)
> (** I need to overcome the above hurdle before evaluating: **)
> logisticnlm = NonlinearModelFit[popvalues2D, Logistic[t, r], {r}, t]
> (** I want to use the initial point Subscript[P, 1] and end point \
> Subscript[P, 49] as "pivot points" and use NonlinearModelFit to get \
> the Best Fit Non-Linear Regression via a "dance a la Levenberg-Marquardt" \
> similar to the dance shown here:
> http://www.numerit.com/samples/nlfit/doc.htm **)
> (** Thank you for your help! **)
The model given to NonlinearModelFit has 2 arguments, but Logistic is only defined for 1 argument. If we use Logistic[t], the estimation works fine.
logisticnlm = NonlinearModelFit[popvalues2D, Logistic[t], {r}, t]
Show[plt1, Plot[logisticnlm[t], {t, 1, 49}]]
You could simplify the model if you wish.
logisticnlm2 = NonlinearModelFit[popvalues2D, Logistic[t] // Simplify, {r}, t]
Show[plt1, Plot[logisticnlm2[t], {t, 1, 49}]]
A look at the difference of the 2 fitted functions shows that they match to within numerical error, though.
Plot[logisticnlm2[t] - logisticnlm[t], {t, 1, 49}]
Darren Glosemeyer
Wolfram Research
- References:
- Dancing "a la Levenberg-Marquardt" to get the best Logistic Model.
- From: Gilmar Rodriguez-pierluissi <peacenova@yahoo.com>
- Dancing "a la Levenberg-Marquardt" to get the best Logistic Model.