Re: [R] integrate dmvtnorm

From: Christos Argyropoulos <argchris_at_hotmail.com>
Date: Thu, 24 Jun 2010 02:43:53 +0300

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.
                                               


        [[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. Received on Thu 24 Jun 2010 - 01:14:05 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:20:34 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