SCM

SCM Repository

[matrix] Diff of /pkg/R/lmer.R
ViewVC logotype

Diff of /pkg/R/lmer.R

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 711, Thu May 5 04:12:44 2005 UTC revision 717, Sun May 8 18:58:05 2005 UTC
# Line 545  Line 545 
545                                              mu = mu,                                              mu = mu,
546                                              wt = weights^2))/2                                              wt = weights^2))/2
547    
548                if (is.null(getOption("laplaceinR")))  ##               if (is.null(getOption("laplaceinR")))
549                {  ##               {
550                    ans <- ans +                    ans <- ans +
551                        .Call("lmer_laplace_devComp", reducedObj,                        .Call("lmer_laplace_devComp", reducedObj,
552                              PACKAGE = "Matrix")                              PACKAGE = "Matrix")
553                }  ##               }
554                else  ##               else
555                {  ##               {
556                    ranefs <- .Call("lmer_ranef", reducedObj, PACKAGE = "Matrix")  ##                   ranefs <- .Call("lmer_ranef", reducedObj, PACKAGE = "Matrix")
557                    ## ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega  ##                   ## ans <- ans + reducedObj@devComp[2]/2 # log-determinant of Omega
558    
559                    Omega <- reducedObj@Omega  ##                   Omega <- reducedObj@Omega
560                    for (i in seq(along = ranefs))  ##                   for (i in seq(along = ranefs))
561                    {  ##                   {
562                        ## contribution for random effects (get it working,  ##                       ## contribution for random effects (get it working,
563                        ## optimize later)  ##                       ## optimize later)
564                        ## symmetrize RE variance  ##                       ## symmetrize RE variance
565                        Omega[[i]] <- Omega[[i]] + t(Omega[[i]])  ##                       Omega[[i]] <- Omega[[i]] + t(Omega[[i]])
566                        diag(Omega[[i]]) <- diag(Omega[[i]]) / 2  ##                       diag(Omega[[i]]) <- diag(Omega[[i]]) / 2
567    
568                        ## want log of `const det(Omega) exp(-1/2 b'  ##                       ## want log of `const det(Omega) exp(-1/2 b'
569                        ## Omega b )` i.e., const + log det(Omega) - .5  ##                       ## Omega b )` i.e., const + log det(Omega) - .5
570                        ## * (b' Omega b)  ##                       ## * (b' Omega b)
571    
572                        ## FIXME: need to adjust for sigma^2 for appropriate  ##                       ## FIXME: need to adjust for sigma^2 for appropriate
573                        ## models (easy).  These are all the b'Omega b,  ##                       ## models (easy).  These are all the b'Omega b,
574                        ## summed as they eventually need to be.  Think of  ##                       ## summed as they eventually need to be.  Think of
575                        ## this as sum(rowSums((ranefs[[i]] %*% Omega[[i]])  ##                       ## this as sum(rowSums((ranefs[[i]] %*% Omega[[i]])
576                        ## * ranefs[[i]]))  ##                       ## * ranefs[[i]]))
577    
578                        ranef.loglik.det <- nrow(ranefs[[i]]) *  ##                       ranef.loglik.det <- nrow(ranefs[[i]]) *
579                            determinant(Omega[[i]], logarithm = TRUE)$modulus/2  ##                           determinant(Omega[[i]], logarithm = TRUE)$modulus/2
580    
581  #                      print(ranef.loglik.det)  ## #                      print(ranef.loglik.det)
582    
583                        ranef.loglik.re <-  ##                       ranef.loglik.re <-
584                            -sum((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]])/2  ##                           -sum((ranefs[[i]] %*% Omega[[i]]) * ranefs[[i]])/2
585    
586  #                      print(ranef.loglik.re)  ## #                      print(ranef.loglik.re)
587    
588                        ranef.loglik <- ranef.loglik.det + ranef.loglik.re  ##                       ranef.loglik <- ranef.loglik.det + ranef.loglik.re
589    
590                        ## Jacobian adjustment  ##                       ## Jacobian adjustment
591                        log.jacobian <-  ##                       log.jacobian <-
592                            sum(log(abs(apply(reducedObj@bVar[[i]],  ##                           sum(log(abs(apply(reducedObj@bVar[[i]],
593                                              3,  ##                                             3,
594    
595                                              ## next line depends on  ##                                             ## next line depends on
596                                              ## whether bVars are variances  ##                                             ## whether bVars are variances
597                                              ## or Cholesly factors  ##                                             ## or Cholesly factors
598    
599                                              ## function(x) sum(diag(x)))  ##                                             ## function(x) sum(diag(x)))
600  ## Was this a bug?                          function(x) sum(diag( La.chol( x ) )))  ## ## Was this a bug?                          function(x) sum(diag( La.chol( x ) )))
601                                              function(x) prod(diag( La.chol( x ) )))  ##                                             function(x) prod(diag( La.chol( x ) )))
602                                        )))  ##                                       )))
603    
604  #                      print(log.jacobian)  ## #                      print(log.jacobian)
605    
606    
607                        ## the constant terms from the r.e. and the final  ##                       ## the constant terms from the r.e. and the final
608                        ## Laplacian integral cancel out both being:  ##                       ## Laplacian integral cancel out both being:
609                        ## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)  ##                       ## ranef.loglik.constant <- 0.5 * length(ranefs[[i]]) * log(2 * base::pi)
610    
611                        ans <- ans + ranef.loglik + log.jacobian  ##                       ans <- ans + ranef.loglik + log.jacobian
612                    }  ##                   }
613                }  ##               }
614                ## ans is (up to some constant) log of the Laplacian                ## ans is (up to some constant) log of the Laplacian
615                ## approximation of the likelihood. Return it's negative                ## approximation of the likelihood. Return it's negative
616                ## to be minimized                ## to be minimized
# Line 626  Line 626 
626    
627            if (method == "Laplace")            if (method == "Laplace")
628            {            {
629  ### Rprof() # trying to figure out if C-ifying bhat is worthwhile  ###Rprof("/tmp/Laplace-profile.out") # trying to figure out if C-ifying bhat is worthwhile
630                ## no analytic gradients or hessians                ## no analytic gradients or hessians
631                optimRes <-                optimRes <-
632                    optim(fn = devLaplace,                    optim(fn = devLaplace,
# Line 691  Line 691 
691            ## have the 'correct' random effects in reducedObj.            ## have the 'correct' random effects in reducedObj.
692    
693            loglik <- devLaplace(optpars)            loglik <- devLaplace(optpars)
694            print(loglik)            ##print(loglik)
695            ff <- optpars[1:(responseIndex-1)]            ff <- optpars[1:(responseIndex-1)]
696            names(ff) <- names(fixef(obj))            names(ff) <- names(fixef(obj))
697    

Legend:
Removed from v.711  
changed lines
  Added in v.717

root@r-forge.r-project.org
ViewVC Help
Powered by ViewVC 1.0.0  
Thanks to:
Vienna University of Economics and Business Powered By FusionForge