Hi Hannah: You have a few things going on, but the bottom line is that numer and denom are both double integrals.
On Thu, Feb 17, 2011 at 1:06 PM, li li <hannah....@gmail.com> wrote: > Hi all, > I have some some problem with regard to finding the integral of a > function containing an indicator function. > please see the code below: > > > func1 <- function(x, mu){ > (mu^2)*dnorm(x, mean = mu, sd = 1)*dgamma(x, shape=2)} > curve(func1(x, 10), 0, 20) curve(func1(x, 5), 0, 20) both yield reasonable plots, so this is a function of one variable when mu is supplied. If you wanted to plot this as a function of mu, you could get a single curve by fixing x or integrating over x, which is what m1star appears to be doing. m1star <- function(x){ > integrate(func1, lower = 0, upper = Inf,x)$val} > u <- seq(0, 20, 0.05) v <- sapply(u, m1star) plot(u, v, type = 'l') yields a plot of m1star, which appears to marginalize func1 to make it a function of mu, from what I can tell. T <- function(x){ > 0.3*dnorm(x)/(0.3*dnorm(x)+0.7*m1star(x))} > plot(u, sapply(u, T), type = 'l') yields a plot of T, but having to use sapply() indicates that T also marginalizes a 2D function. > func2 <- function(x,c){(T(x) <=c)*0.3*dnorm(x)} > func3 <- function(x,c){(T(x) <= c)*(0.3*dnorm(x)+0.7*m1star(x))} > func2 uses T, which in turn uses m1star, so func2 is a marginalization of a 2D function whose support is on T(x) <= c. To integrate func2, you have to do the integration in m1star first, so basically you have a double integral to evaluate in numer. The same problem occurs in func3, since m1star() is a part of both T() and the convex combination. Therefore, denom also requires double integration. > numer <- function(c){ > integrate(func2, -Inf, Inf, c)$val} > > denom <- function(c){ > integrate(func3, lower, Inf,c)$val} > > Look into cubature or Rcuba for two packages that are capable of performing numerical double integration. An alternative solution is to use approxfun() to create a function object from evaluations of m1star and T, and then use the approxfun()s as the functions to be integrated in numer and denom. If you decide to go the approxfun route, make sure to make enough evaluations to reasonably capture the curvature in both m1star and T. HTH, Dennis > > > The error message is as below : > > > numer(0.5) > Error in integrate(func1, lower = 0, upper = Inf, x) : > the integral is probably divergent > > [[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. > [[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.