> Dear all, > > How can I fit a power model in R. > > Thanks in advance, > > miltinho > > > > para armazenamento! > > [[alternative HTML version deleted]]
Try the approach below with the lognormal distribution for the data. The data is real. Rubén y<-c(2841,1151,1579,1491,1306,2294,1781,1472,2709,1845,1882,2113,3683,2852,794,1965,3417,4478,3523, 2414,2769,5216,3590,4057,2770,3171,3208,6265,4411,2319,3044,1771,3472,4394,3169,4225,4934,3916,4120, 5168,4461,5980,5670,5806,5479,3644,5529,5011,7826,3964,6942,4285,4301,4696,4288,4921,7138, 9881,5025,5181,7110,5413,8279,7803,4652,5114,3127,5798,4048,7459,4741,5798,6738,3562,6235,8015, 5140,5852,4634,5006,5937,5014,9608,2134,7649,7649,6634,5028,5943,4732,7782,7035,8147,6583,5085, 5743,10932,6935,5716,5716,5839,6356,10398,9709,11937,5955,9385,17764,15560,13834,20396,20961) x<-c(8.0,8.5,8.5,9.0,9.0,9.0,10.0,10.8,11.0,11.0,11.2,11.5,11.5,11.5,11.6,12.0,12.0,12.0,12.0,12.5,12.5,12.6, 12.8,12.8,12.9,13.0,13.0,13.5,13.5,13.5,13.7,14.0,14.0,14.0,14.0,14.2,14.5,14.5,14.5,14.8,15.0,15.0,15.5,15.5, 15.5,15.5,15.9,16.0,16.0,16.0,16.0,16.0,16.0,16.5,16.5,16.5,16.5,16.5,16.5,16.5,16.5,16.5,16.5,16.5,16.6,17.0, 17.0,17.0,17.0,17.0,17.0,17.0,17.0,17.5,17.5,17.5,17.5,17.5,17.5,17.5,17.7,18.0,18.0,18.0,18.0,18.0,18.0,18.0, 18.0,18.5,18.5,19.0,19.0,19.0,19.0,19.0,19.5,19.5,19.5,19.5,19.5,19.5,20.0,20.5,21.0,22.0,24.5,26.5,26.5,26.8, 28.0,30.0) f1=2 #starting value for parameter f0=20 #starting value for parameter y.mod<-f0*(x)^f1 plot(x,y) lines(x,y.mod) fn<-function(p){ y.mod<-p[1]*x^p[2]; squdiff<-(log(y)-log(y.mod))^2; lik=(length(y)/2)*log(sum(squdiff)/length(y)) #lognormal } y.lik<-nlm(fn,p=c(20,2),hessian=TRUE) y.lik f0<-y.lik$estimate[1] f1<-y.lik$estimate[2] covmat<-solve(y.lik$hessian) #estimated covariance matrix covmat sef0<-sqrt(covmat[1,1]) #standard error sef1<-sqrt(covmat[2,2]) #standard error covf0f1<-covmat[1,2] sef0 sef1 y.fit<-f0*x^f1 plot(x,y,pch=19,xlab="",ylab="") lines(x,y.fit,lwd=par(3)) text(x=10,y=20000,expression(paste("y=",19.9*x^1.98))) ______________________________________________ R-help@r-project.org mailing list 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.