[R] Problems with adapt
Cristian Angelo Guevara
crguevar at mit.edu
Mon Mar 23 15:30:43 CET 2009
Hi:
I'm trying to estimate a model which involves the estimation of double
integrals, so I'm using adapt procedure. I need to calculate the
integrals trough my 2000 size database, so I do it using a loop. My
code correctly estimates the integral for the first row, but for the
second R crashes. I tried changing the order of the data but the
result is the same so I guess there is something wrong with my code of
the loop or, more probably, on the use of adapt. I'm new to R so I
would appreciate any comment on my code below.
Thanks
Angelo
----------------------------------------------------------------------------------------------------------
DATA<-read.table("Data.csv",header=TRUE,sep=",")
library(adapt)
mnl.lik<-function(theta,y){
th1<-theta[1]
th2<-theta[2]
tha<-theta[3]
thb<-theta[4]
thc<-theta[5]
thp<-theta[6]
thmu<-theta[7]
alfz1<-theta[8]
alfz2<-theta[9]
alfc<-theta[10]
alf<-theta[11]
r<- 1
s<- 1
n<-2000
lik<-numeric(n)
int<-numeric(n)
v<- numeric(2)
for (i in 1:n) { #Beggin Loop
lstarpre<- function(v){ #This is the fuction to be integrated
e1<-y$p1[i]-alfz1*y$z_a1[i] - alfz2*y$z_b1[i] -alfc*y$c1[i] -alf
e2<-y$p2[i]-alfz1*y$z_a2[i] - alfz2*y$z_b2[i] -alfc*y$c2[i] -alf
e3<-y$p3[i]-alfz1*y$z_a3[i] - alfz2*y$z_b3[i] -alfc*y$c3[i] -alf
U1<- (th1 +tha*y$a1[i] +thb*y$b1[i] +thc*y$c1[i] +thp*y$p1[i]
+thmu*(e1[i]+v[1]))
U2<- (th2 +tha*y$a2[i] +thb*y$b2[i] +thc*y$c2[i] +thp*y$p2[i]
+thmu*(e2[i]+v[2]))
U3<- ( +tha*y$a3[i] +thb*y$b3[i] +thc*y$c3[i] +thp*y$p3[i] +thmu*(e3[i]))
Pni<- (y$Ch1[i]*exp(U1) + y$Ch2[i]*exp(U2)+
y$Ch3[i]*exp(U3))/(exp(U1)+exp(U2)+exp(U3))
fe<- exp((-e1[i]*e1[i]-e2[i]*e2[i]-e3[i]*e3[i])/(2*r))
fv<- exp((-v[1]*v[1]-v[2]*v[2])/(2*s))
prelik<- Pni*fe*fv
return(prelik)
}
int[i]<- adapt(2,lo = c(-4,-4), up = c(4,4), functn= lstarpre, eps=0.01)$value
print(int[i])
lik[i]<- log(max(1E-20,int))
print(lik[i])
} #End Loop
logl<-sum(lik)
return(-logl)
}
p<-optim(c(0,0,0,0,0,0,0,0,0,0,0),mnl.lik,y=DATA,method="BFGS",hessian=T,
control = list(maxit=6000, temp=2000, trace=4))
More information about the R-help
mailing list