[R] Power model in R

[Ruben Roa Ureta]

> 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)))