[R] integrate dmvtnorm

[Christos Argyropoulos]

Carrie, 

According to the help file for the function "integrate",  the integrand 'f': 

     must accept a vector of inputs and produce a vector of
     function evaluations at those points.  The ‘Vectorize’ function
     may be helpful to convert ‘f’ to this form.

The integrate function wraps the QUADPACK library of numerical integrators 
(http://www.netlib.org/quadpack/readme) most (if not all) of which are 
interpolatory quadrature rules (preceded by a variable transformation for 
infinite domains of integration).  Therefore they approximate your integral by
a finite weighted sum using a set of weights and function evaluations at an 
equally size set of nodes. Conceptually (although I doubt that it actually works
this way inside QUADPACK), one can simply form the inner product of the 
vector of weights and the vector of function values at the nodes and use this
as an estimate of the value of the integral. For an adaptive rule one repeats 
this procedure by increasing the number of knots until the desired accuracy 
has been achieved. From this simple description you can see why the integrand
has to be a "vectorized" function (the integrator will ask 'f' for a vector of function
evaluations at the set of nodes). This can be achieved either explicitly or implicitly 
(by calling the sapply inside the function).


Christos  


Date: Wed, 23 Jun 2010 20:22:15 -0400
Subject: Re: [R] integrate dmvtnorm
From: carrieandstat at gmail.com
To: argchris at hotmail.com
CC: r-help at r-project.org

Thanks! 

Both suggestions are very helpful.

One more question:

Can I use Vectorize to solve double integration question ?



Now that 

f=function(x, y) {dnorm(y, mean= 0.75/x)*dnorm(x, 
mean=0.6,

 sd=0.15)}



And I want to integrate x first,then y. Ravi used sapply, which is good,
 but it seems to be that Vectorize is easier.



Thanks for help. I appreciated



Carrie

2010/6/23 Christos Argyropoulos <argchris at hotmail.com>






No something else is going on here ....

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

> f(1)
[1] 0.01194131

> x<-seq(-2,2,.15)
> f(x)
Error in dmvnorm(c(0.6, 0.8), mean = c(0.75, 0.75/x)) : 
  mean and sigma have non-conforming size


But ...

> sapply(x,f)
 [1] 1.205791e-66 2.377822e-59 1.712003e-52 4.488794e-46 4.269526e-40
 [6] 1.464321e-34 1.793031e-29 7.702766e-25 1.122712e-20 5.165600e-17
[11] 6.242351e-14 1.074366e-11 8.904914e-12 2.165575e-59 2.892453e-13

[16] 2.446326e-03 9.655456e-02 3.377855e-01 3.230318e-01 1.040144e-01
[21] 1.194131e-02 4.984067e-04 7.620137e-06 4.281072e-08 8.849889e-11
[26] 6.735400e-14 1.887638e-17



suggesting the solution:

vf<-Vectorize(f)

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


Christos




> Date: Wed, 23 Jun 2010 19:05:53 -0400

> From: carrieandstat at gmail.com
> To: R-help at r-project.org
> 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.

 		 	   		  
Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign up now.

 		 	   		  
_________________________________________________________________
Hotmail: Powerful Free email with security by Microsoft.