integrate dmvtnorm

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

integrate dmvtnorm

Carrie Li
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]]

______________________________________________
[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.
Reply | Threaded
Open this post in threaded view
|

Re: integrate dmvtnorm

Ravi Varadhan
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: [hidden email] [mailto:[hidden email]] 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]]

______________________________________________
[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.
Reply | Threaded
Open this post in threaded view
|

Re: integrate dmvtnorm

Christos Argyropoulos
In reply to this post by Carrie Li

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


        [[alternative HTML version deleted]]

______________________________________________
[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.
Reply | Threaded
Open this post in threaded view
|

Re: integrate dmvtnorm

Carrie Li
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 <[hidden email]>

>  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: [hidden email]
> > To: [hidden email]
> > 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]]
> >
> > ______________________________________________
> > [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.
>
> ------------------------------
> Your E-mail and More On-the-Go. Get Windows Live Hotmail Free. Sign up
> now. <https://signup.live.com/signup.aspx?id=60969>
>

        [[alternative HTML version deleted]]

______________________________________________
[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.
Reply | Threaded
Open this post in threaded view
|

Re: 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: [hidden email]


----- Original Message -----
From: Carrie Li <[hidden email]>
Date: Wednesday, June 23, 2010 8:23 pm
Subject: Re: [R] integrate dmvtnorm
To: Christos Argyropoulos <[hidden email]>
Cc: [hidden email]


> 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 <[hidden email]>
>  
>  >  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: [hidden email]
>  > > To: [hidden email]
>  > > 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]]
>  > >
>  > > ______________________________________________
>  > > [hidden email] 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]]
>  
>  ______________________________________________
>  [hidden email] mailing list
>  
>  PLEASE do read the posting guide
>  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.
Reply | Threaded
Open this post in threaded view
|

Re: integrate dmvtnorm

Christos Argyropoulos
In reply to this post by Carrie Li

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: [hidden email]
To: [hidden email]
CC: [hidden email]

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






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: [hidden email]
> To: [hidden email]
> 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]]

>
> ______________________________________________
> [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.

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

     
_________________________________________________________________
Hotmail: Powerful Free email with security by Microsoft.

______________________________________________
[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.