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.hlx@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]]