[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