[R] integrate dmvtnorm

[Ravi Varadhan]

Double integral:

fvec = function(x, y) sapply(x, function(z, y) dnorm(y,  mean=0.75/z) * dnorm(z, mean=0.6, sd=0.15), y=y)  

gvec = function(x) sapply(x, function(y) integrate(f, lower=-Inf, upper=Inf, y=y, subdivisions=10000, rel.tol=1.e-08)$val)

xx <- seq(-5, 5, length=1000)

plot(xx, gvec(xx), type="l")

> integrate(gvec, lower=-Inf, upper=Inf, subdivisions=10000, rel.tol=1.e-06)
0.999914 with absolute error < 5.7e-07

> integrate(gvec, lower=0, upper=Inf, subdivisions=10000, rel.tol=1.e-06)
0.8954554 with absolute error < 9.4e-07
> 

Ravi.
____________________________________________________________________

Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology
School of Medicine
Johns Hopkins University

Ph. (410) 502-2619
email: rvaradhan at jhmi.edu


----- Original Message -----
From: Carrie Li <carrieandstat at gmail.com>
Date: Wednesday, June 23, 2010 8:23 pm
Subject: Re: [R] integrate dmvtnorm
To: Christos Argyropoulos <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
>  > > 
>  > > PLEASE do read the posting guide
>  > 
>  > > and provide commented, minimal, self-contained, reproducible code.
>  >
>  > ------------------------------
>  > Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign 
> up
>  > now. <>
>  >
>  
>  	[[alternative HTML version deleted]]
>  
>  ______________________________________________
>  R-help at r-project.org mailing list
>  
>  PLEASE do read the posting guide 
>  and provide commented, minimal, self-contained, reproducible code.