интегрирование по интеграл в R

Question

Я хочу решить следующие в R:

∫ ^Н [π ( т) ∫ _т^Н A ( x) dx] dt

Где π (t) - предшествующее, а A (x) - это функция A, определенная ниже.

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

С π (т), А (х)> 0, expected_loss (0,5) должно быть меньше, чем expected_loss (1). Но это не то, что я получаю:

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

Я не уверен, что я делаю неправильно.

Answer 1

В вашем integrand, lower = t не является векторизованным, поэтому вызов для интеграции не делает то, что вы ожидали *. Vectorising более t исправляет эту проблему,

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 

*: пристальный взгляд на integrate показал, что переход векторов понизить и/или верхней молча допускается, но только первое значение было принято во внимание. При интегрировании на более широком интервале квадратурная схема выбрала первую точку дальше от начала координат, что привело к неинтуитивному уменьшению, которое вы наблюдали.

После сообщения об этом поведении r-devel, this user-error will now be caught by integrate благодаря Martin Maechler (R-devel).

Answer 2

В данном конкретном случае вам не нужно Vectorize, так как интеграл от dbeta уже реализован в R до pbeta. Попробуйте следующее:

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