Tích hợp trên một không thể thiếu trong R

Question

Tôi muốn giải quyết sau trong R:

∫ ^H [π ( t) ∫ _t^H A ( x) dx] dt

Trường hợp π (t) là trước và A (x) là hàm A được xác định bên dưới.

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

Kể từ π (t), A (x)> 0, expected_loss (0,5) nên được ít hơn expected_loss (1). Nhưng đây không phải là những gì tôi nhận được:

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

Tôi không chắc mình đang làm gì sai.

Answer 1

Trong số integrand, lower = t của bạn không được vectơ hóa, vì vậy cuộc gọi tích hợp không thực hiện những gì bạn mong đợi *. Vectorising qua t sửa chữa vấn đề này,

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 

*: một cái nhìn sâu hơn về integrate tiết lộ rằng qua vectơ để giảm và/hoặc trên đã âm thầm cho phép, nhưng chỉ có giá trị đầu tiên được đưa vào tính toán. Khi tích hợp trong khoảng thời gian rộng hơn, sơ đồ phương trình đã chọn điểm đầu tiên xa hơn điểm gốc, dẫn đến giảm không trực quan mà bạn quan sát được.

Sau khi báo cáo hành vi này cho r-devel, this user-error will now be caught by integrate nhờ Martin Maechler (R-devel).

Answer 2

Trong trường hợp cụ thể này, bạn không cần phải Vectorize vì tích phân của dbeta đã được triển khai trong R qua pbeta. Hãy thử điều này:

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