[R] Generation of Multariate Gumbel Random Numbers VGAM or GUMBEL

2009-09-27 Thread Cristian Angelo Guevara
Dear R users:

I'm trying to generate multivariate random numbers where some
variables are correlated. I tried VGAM and GUMBEL, but I really don't
understand from the documentation, how to move from the single random
to the multivariate. Does anybody have any insight about about this?

Thanks in advance
Angelo

__
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.


[R] Problems with adapt

2009-03-23 Thread Cristian Angelo Guevara
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))

__
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.


[R] integration within maximum likelihood

2008-12-04 Thread Cristian Angelo Guevara
Hi:

I'm trying to estimate a latent variable model in mnl discrete choice
framework using R. I need to do first a uni dimensional integral
within each observation (row) in the database and then sum over
observations. I'm stacked in the point shown below. Apparently I have
a dimensionality problem in the definition of the integral. Maybe it
does not identify that what I need is only one row at a time. I'm no
sure.

I would appreciate any kind of help. Below is my code and the error report.

Thanks
Angelo


> setwd("C:/Users/Angelo/Documents/PhD/Montecarlo TVI")
> DATA<-read.table("Data.csv",header=TRUE,sep=",")
> 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]
+ alfz<-theta[8]
+ alfp<-theta[9]
+ mu1<- alfz*y$z_a1 + alfp*y$p1
+ mu2<- alfz*y$z_a2 + alfp*y$p2
+ mu3<- alfz*y$z_a3 + alfp*y$p3
+ U1<- th1 +tha*y$a1 +thb*y$b1 +thc*y$c1  +thp*y$p1 +thmu*mu1
+ U2<- th2 +tha*y$a2 +thb*y$b2 +thc*y$c2  +thp*y$p2 +thmu*mu2
+ U3<- +tha*y$a3 +thb*y$b3 +thc*y$c3  +thp*y$p3 +thmu*mu3
+ Usum<- (exp(U1)+exp(U2)+exp(U3))
+ arg1<-  function(mu1) {(exp(U1)/Usum)*dnorm(mu1- alfz*y$z_a1 -alfp*y$p1)}
+ arg2<-  function(mu2) {(exp(U2)/Usum)*dnorm(mu2- alfz*y$z_a2 -alfp*y$p2)}
+ arg3<-  function(mu3) {(exp(U3)/Usum)*dnorm(mu3- alfz*y$z_a3 -alfp*y$p3)}
+ int1<- integrate(arg1,-Inf,Inf)
+ int2<- integrate(arg2,-Inf,Inf)
+ int3<- integrate(arg3,-Inf,Inf)
+ logl<-sum(y$Ch1*int1$value+ y$Ch2*int2$value + y$Ch3*int3$value)
+ return(-logl)
+ }
>
> p<-optim(c(0,0,0,0,0,0,0,0,0),mnl.lik,y=DATA,method="BFGS",hessian=T)
Error in integrate(arg1, -Inf, Inf) :
  evaluation of function gave a result of wrong length
Calls: optim ->  -> fn -> integrate -> .External
In addition: Warning message:
In mu1 - alfz * y$z_a1 :
  longer object length is not a multiple of shorter object length
Execution halted

__
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.