Integrate over an integral in R -
i want solve following in r:
∫0h [π(t) ∫th a(x) dx ] dt
where π(t) prior , a(x) function defined below.
prior <- function(t) dbeta(t, 1, 24) <- 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 less expected_loss(1). not get:
> expected_loss(.5) [1] 0.2380371 > expected_loss(1) [1] 0.0625
i'm not sure i'm doing wrong.
in integrand
, lower = t
not vectorised, call integrate not doing expected*. vectorising on t
fixes 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
*: closer @ integrate
revealed passing vectors lower and/or upper silently allowed, first value taken account. when integrating on wider interval quadrature scheme picked first point further origin, resulting in unintuitive decrease observed.
after reporting behaviour r-devel, this user-error caught integrate martin maechler (r-devel).
Comments
Post a Comment