Announcement

Collapse
No announcement yet.
X
  • Filter
  • Time
  • Show
Clear All
new posts

  • double integral in Mata

    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!







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

    Comment


    • #3
      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!

      Comment


      • #4
        If u1 and u2 are jointly normal then perhaps check out the multivariate normal integral functions at
        Code:
        help mf_mvnormal

        Comment


        • #5
          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.

          Comment


          • #6
            Originally posted by Tianji Cai View Post
            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.

            Comment

            Working...
            X