I want to solve the following in R:
∫0H [π(t) ∫tHA(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.
In short, you may use R to find out a numerical answer to an n-fold integral. −5 . For more information about the function integrate, type help(integrate) in R.
Since integration is the inverse of differentiation, many differentiation rules lead to corresponding integration rules. Consider, for example, the chain rule. The formula forms the basis for a method of integration called the substitution method.
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).
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
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With