Gaussian quadrature for bivariate normal distribution

[Niroshan]

Dear R Users

I am maximizing  a likelihood function it has two two correlated random effects which follows bivariate normal distribution. To get the marginal distribution I want to integrate out with respect to these two correlated random effects. Does any body know how can I implement gaussian quadrature to approximate this integral.

Thanks for taking time

[Michael Weylandt]

What distribution's log-likelihood are you using? If sigma is supposed
to be square-rooted, you may have to put a constraint or use of abs()
might suffice -- though, admittedly I'm not sure what that will do to
convergence behavior -- but it's hard to help without seeing the
function at hand.

Michael

On Wed, Mar 14, 2012 at 6:57 PM, Niroshan <[hidden email]> wrote:

> Dear R Users
>
>
> I am maximizing a user defined log likelihood function. It includes variance
> parameter (sigma). I used R function optim with BFGS maximization method.
> However, it stops before the solution saying “sqrt(sigma): NaNs produced”
> Could anybody know a proper transformation for sigma which can be passed in
> the function?  For the correlation parameter I used Fishers’ transformation
> so it didn’t give any problems. I know we can use ‘L-BFGS-B’ method instead
> of ‘BFGS’. However, it produces some non sensible estimates.
>
> Thanks for taking time
>
>
> --
> View this message in context: http://r.789695.n4.nabble.com/Maximization-problem-in-the-optim-function-tp4473408p4473408.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> [hidden email] 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.

______________________________________________
[hidden email] 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.

[Niroshan]

This is the part of the function in my likelihood function




prob<-function(t1,t2,m1,m2,s_r,s_p,s_e,cor_m){

if ((t1==0) & (t2==0))

log_lik<-log(pmvnorm(lower=rep(-Inf, 2), upper=c(m1/sqrt(s_r+s_p+s_e+1),m2/sqrt(s_r+s_p+s_e+1)),
 mean=rep(0,2), corr=cor_m))

else if ((t1==0) & (t2==1))

log_lik<-log(pmvnorm(lower=c(-Inf,m2/sqrt(s_r+s_p+s_e+1)), upper=c(m1/sqrt(s_r+s_p+s_e+1),Inf),
 mean=rep(0,2), corr=cor_m))

else if ((t1==1) & (t2==0))

log_lik<-log(pmvnorm(lower=c(m1/sqrt(s_r+s_p+s_e+1),-Inf), upper=c(Inf,m2/sqrt(s_r+s_p+s_e+1)),
 mean=rep(0,2), corr=cor_m))

else if ((t1==1) & (t2==1))

log_lik<-log(pmvnorm(lower=c(m1/sqrt(s_r+s_p+s_e+1),m2/sqrt(s_r+s_p+s_e+1)), upper=c(Inf,Inf),
 mean=rep(0,2), corr=cor_m))

return(log_lik)
}