Re: [R] integrate dmvtnorm

From: Ravi Varadhan <rvaradhan_at_jhmi.edu>
Date: Wed, 23 Jun 2010 19:26:52 -0400

The main problem is that your function is not vectorized. Here is one solution:

> require(mvtnorm)

> f=function(x){ sapply(x, function(y) {dmvnorm(c(0.6, 0.8), mean=c(0.75,
0.75/y))*dnorm(y, mean=0.6, sd=0.15)}) }

> integrate(f, lower=-Inf, upper=Inf)
0.1314427 with absolute error < 4e-05
>

Hope this helps,
Ravi.

-----Original Message-----
From: r-help-bounces_at_r-project.org [mailto:r-help-bounces_at_r-project.org] On Behalf Of Carrie Li
Sent: Wednesday, June 23, 2010 7:06 PM
To: r-help
Subject: [R] integrate dmvtnorm

Hello, everyone,

I have a question about integration of product of two densities. Here is the sample code; however the mean of first density is a function of another random variable, which is to be integrated.

##
f=function(x) {dmvnorm(c(0.6, 0.8), mean=c(0.75, 0.75/x))*dnorm(x, mean=0.6, sd=0.15)}
integrate(f, lower=-Inf, upper=Inf)

## error message
Error in dmvnorm(c(0.6, 0.8), mean = c(0.75, 0.75/x)) :   mean and sigma have non-conforming size

I think it's because the mean in dmvnorm is a function of x....

is there any package or function to handle this question ?

Thanks for any help!

Carrie

        [[alternative HTML version deleted]]



R-help_at_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-help_at_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. Received on Wed 23 Jun 2010 - 23:28:45 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Thu 24 Jun 2010 - 01:10:36 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.

list of date sections of archive