Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Integrate over an integral in R

Tags:

r

integral

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.

like image 717
LJB Avatar asked Mar 29 '16 05:03

LJB


People also ask

Can R solve integrals?

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.

Does chain rule work with integration?

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.


2 Answers

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).

like image 163
baptiste Avatar answered Nov 02 '22 23:11

baptiste


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
like image 30
nicola Avatar answered Nov 03 '22 00:11

nicola