double integral in Mata

[Tianji Cai]

Hey all,

I am working on to implement a random effect zip model in Mata. The model includes two random effects: one in the logit part, and another is in poisson part. Just wondering if anyone could give me an example how to do a double integration in Mata.

For example, my cluster level likelihood is as follows.



real colvector hiu1u2(real colvector u1, real colvector u2, real colvector Ai, real colvector Bi, real colvector zi, real colvector yi, real scalar s1, real scalar s2, real scalar t) {
n1=rows(u1)
n2=rows(u2)
r=J(n1*n2,1,0)
for(i=1; i<=n1;++i){
eta=Ai+s1*u1[i]
for(j=1;j<=n2;++j){
zeta=Bi+s2*u2[j]
lambda=exp(zeta)
pi=1/(1+exp(-eta))
p0=zi:*log(exp(eta)+exp(-lambda))+(1:-zi)*(-lambda+yi:*zeta-lngamma(yi+1))+log(1+exp(eta))
r[(i-1)*n1+j]=p0-.5*u1[i]^2-.5*u2[j]^2
}
return(r)
}


Just dont know how to conduct the double integration over u1 & u2. I tried the Quadrature but had problem to give upper and lower bounds for u1 and u2 simultaneously.

Any suggestion will be highly appreciated!

Thank you very much!

[Joseph Coveney]

That resembles code used with SAS PROC NLMIXED . Have you looked into whether you can fit this model with menl?

[Tianji Cai]

Yes, SAS can do it. I want to modify the random zip model by adding some features. So, I really need to implement it in Mata. Thanks!

[John Mullahy]

If u1 and u2 are jointly normal then perhaps check out the multivariate normal integral functions at

help mf_mvnormal

[Tianji Cai]

Thank you very much! I figured out how to use Quadrature to get it done. I will post the code if anyone else is interested once the testing is done.

[Joseph Coveney]

I figured out how to use Quadrature to get it done. I will post the code if anyone else is interested once the testing is done.

Yes, please do. I am interested.

[Ahmed Sadek Yousuf]

Thank you very much! I figured out how to use Quadrature to get it done. I will post the code if anyone else is interested once the testing is done.

Hi ! I am interested in the code. Can you share ?