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: carrieands...@gmail.com To: argch...@hotmail.com CC: r-help@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 <argch...@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: carrieands...@gmail.com > To: R-help@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@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. ______________________________________________ R-help@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.
