Hi Dennis,
Thank you for your kind reply. Yes, essentially, we take integration
twice. However, I still have a few questions:
First, if we consider doing a double integration at the end, since the
first integration in within the indicator function,
it seems to be difficult.
Second, m1star is a univariate function, say, of x, since I already take
the integration. Then T is an univariate
function of x and func2 is a function of x and c. It seems to make sense
if I integrate func2 with respect to x.
Can you give more hint in terms of doing a double integration? Thank you
again.
Hannah
2011/2/17 Dennis Murphy <djmu...@gmail.com>
> 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<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.
