[R] nls in r
ProfJCNash
profjcnash at gmail.com
Wed Aug 12 20:25:47 CEST 2015
With package nlmrt, I get a solution, but the Jacobian is essentially
singular, so the model may not be appropriate. You'll need to read the
documentation to learn how to interpret the Jacobian singular values. Or
Chapter 6 of my book "Nonlinear parameter optimization with R tools."
Here's the script. Note nlxb is a little different from nls() and
requires a data frame for the data.
library(stats)
x=c(30:110)
y=c(0.000760289, 0.000800320, 0.000830345, 0.000840353, 0.000860370,
0.000910414, 0.000990490, 0.001090594, 0.001200721, 0.001350912,
0.001531172, 0.001751533, 0.001961923, 0.002192402, 0.002463031,
0.002793899, 0.003185067, 0.003636604, 0.004148594, 0.004721127,
0.005394524, 0.006168989, 0.007014544, 0.007870894, 0.008758242,
0.009656474, 0.010565620, 0.011485709, 0.012396520, 0.013308162,
0.014271353, 0.015326859, 0.016525802, 0.017889059, 0.019447890,
0.021223636, 0.023145810, 0.025174229, 0.027391752, 0.029645106,
0.032337259, 0.035347424, 0.039375125, 0.043575783, 0.048003973,
0.052926206, 0.058307309, 0.064581189, 0.071947231, 0.080494476,
0.089891885, 0.100671526, 0.111971207, 0.124237571, 0.137975539,
0.153629149, 0.171239194, 0.190712664, 0.212079979, 0.235026373,
0.259457493, 0.282867017, 0.307830359, 0.334773680, 0.364001719,
0.395742526, 0.425596389, 0.458391314, 0.494558651, 0.534657357,
0.579354317, 0.616075034, 0.656680256, 0.701804548, 0.752133146,
0.808558032, 0.872226001, 0.944664487, 1.027837007, 1.124484096,
1.238426232)
vdata<- data.frame(x=x, y=y)
a=0.0189
b=0.14328
delta=0.0005
# fit = nls(y ~
a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+(a*(exp(b*(x+0.5))-1))))^(0.5),
# start=list(a=a,b=b, delta=delta), trace=TRUE)
# predict(fit)
# plot(x,y,col="red", xlab="Usia",ylab=expression(paste(mu)))
# lines(x,predict(fit), col="blue")
# legend("topleft",
# c(expression(paste(mu)),"Fit"),col=c("red","blue"),lty=1:1)
library(nlmrt)
vformula <- y ~
a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+a*(exp(b*(x+0.5))-1)))^(0.5)
fitjn = nlxb(formula= vformula, start=list(a=a,b=b,delta=delta),
trace=TRUE, data=vdata)
fitjn
Best, JN
On 15-08-12 06:37 AM, vidya wrote:
> I get this error
> Error in numericDeriv(form[[3L]], names(ind), env) :
> Missing value or an infinity produced when evaluating the model
> I was replace the starting value but still get error.
>
> Here is my code:
> library(stats)
> x=c(30:110)
> y=c(0.000760289, 0.000800320, 0.000830345, 0.000840353, 0.000860370,
> 0.000910414, 0.000990490, 0.001090594, 0.001200721, 0.001350912,
> 0.001531172, 0.001751533, 0.001961923, 0.002192402, 0.002463031,
> 0.002793899, 0.003185067, 0.003636604, 0.004148594, 0.004721127,
> 0.005394524, 0.006168989, 0.007014544, 0.007870894, 0.008758242,
> 0.009656474, 0.010565620, 0.011485709, 0.012396520, 0.013308162,
> 0.014271353, 0.015326859, 0.016525802, 0.017889059, 0.019447890,
> 0.021223636, 0.023145810, 0.025174229, 0.027391752, 0.029645106,
> 0.032337259, 0.035347424, 0.039375125, 0.043575783, 0.048003973,
> 0.052926206, 0.058307309, 0.064581189, 0.071947231, 0.080494476,
> 0.089891885, 0.100671526, 0.111971207, 0.124237571, 0.137975539,
> 0.153629149, 0.171239194, 0.190712664, 0.212079979, 0.235026373,
> 0.259457493, 0.282867017, 0.307830359, 0.334773680, 0.364001719,
> 0.395742526, 0.425596389, 0.458391314, 0.494558651, 0.534657357,
> 0.579354317, 0.616075034, 0.656680256, 0.701804548, 0.752133146,
> 0.808558032, 0.872226001, 0.944664487, 1.027837007, 1.124484096,
> 1.238426232)
>
> a=0.0189
> b=0.14328
> delta=0.0005
>
> fit = nls(y ~
> a*(exp(b*(x+0.5)))*((delta*b)/((delta*b)+(a*(exp(b*(x+0.5))-1))))^(0.5),
> start=list(a=a,b=b, delta=delta))
> predict(fit)
> plot(x,y,col="red", xlab="Usia",ylab=expression(paste(mu)))
> lines(x,predict(fit), col="blue")
> legend("topleft",
> c(expression(paste(mu)),"Fit"),col=c("red","blue"),lty=1:1)
>
>
> I really appreciate for the helps. Thank you.
>
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/nls-in-r-tp4711012.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at r-project.org mailing list -- To UNSUBSCRIBE and more, see
> https://stat.ethz.ch/mailman/listinfo/r-help
> PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
> and provide commented, minimal, self-contained, reproducible code.
>
More information about the R-help
mailing list