CODEWITHSUNDEEP
rintegral
LJB
LJB

Reputation: 353

Integrate over an integral in R

I want to solve the following in R:

∫0H [π(t) ∫tH A(x) dx ] dt

Where π(t) is the prior and A(x) is the A function defined below. 

prior <- function(t) dbeta(t, 1, 24)
A     <- function(x) dbeta(x, 1, 4)
expected_loss <- function(H){
  integrand     <- function(t) prior(t) * integrate(A, lower = t, upper = H)$value
  loss          <- integrate(integrand, lower = 0, upper = H)$value
  return(loss)
}

Since π(t), A(x) > 0, expected_loss(.5) should be less than expected_loss(1). But this is not what I get:

> expected_loss(.5)
[1] 0.2380371
> expected_loss(1)
[1] 0.0625

I'm not sure what I'm doing wrong.

Upvotes: 5

Views: 1540

Answers (2)

baptiste
baptiste

Reputation: 77096

In your integrand, lower = t is not vectorised, so the call to integrate is not doing what you expected*. Vectorising over t fixes this issue,

expected_loss <- function(H){
  integrand <- function(t) prior(t) * integrate(A, lower = t, upper = H)$value
  vint <- Vectorize(integrand, "t")
  loss <- integrate(vint, lower = 0, upper = H)$value
  return(loss)
} 

expected_loss(.5)
# [1] 0.7946429
expected_loss(1)
# [1] 0.8571429

*: a closer look at integrate revealed that passing vectors to lower and/or upper was silently allowed, but only the first value was taken into account. When integrating over a wider interval the quadrature scheme picked a first point further from the origin, resulting in the unintuitive decrease that you observed.

After reporting this behaviour to r-devel, this user-error will now be caught by integrate thanks to Martin Maechler (R-devel).

Upvotes: 8

nicola
nicola

Reputation: 24480

In this particular case, you don't need to Vectorize since the integral of dbeta is already implemented in R through pbeta. Try this:

prior <- function(t) dbeta(t, 1, 24)
#define the integral of the A function instead
Aint     <- function(x,H) pbeta(H, 1, 4) - pbeta(x,1,4)
expected_loss <- function(H){
  integrand<-function(x) Aint(x,H)*prior(x)
  loss          <- integrate(integrand, lower = 0, upper = H)$value
  return(loss)
}
expected_loss(.5)
#[1] 0.7946429
expected_loss(1)
#[1] 0.8571429

Upvotes: 6

Related Questions