Integrare su un integrale in R

Question

voglio risolvere il seguente in R:

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

Dove π (t) è il precedente e A (x) è la funzione A definita di seguito.

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

Da π (t), A (x)> 0, expected_loss (.5) dovrebbe essere inferiore a expected_loss (1). Ma questo non è quello che ottengo:

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

Non sono sicuro di quello che sto facendo male.

Answer 1

Nel tuo integrand, lower = t non è vettorizzato, quindi la chiamata per l'integrazione non sta facendo quello che ti aspettavi *. Vectorising oltre t correzioni questo problema,

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 

*: uno sguardo più da vicino integrate ha rivelato che il passaggio vettori per abbassare e/o superiore era silenziosamente consentito, ma solo il primo valore è stato preso in considerazione. Quando si integrava su un intervallo più ampio, lo schema di quadratura preleva un primo punto più lontano dall'origine, determinando la diminuzione non intuitiva osservata.

Dopo aver segnalato questo comportamento a r-devel, this user-error will now be caught by integrate grazie a Martin Maechler (R-devel).

Answer 2

In questo caso particolare, non è necessario Vectorize poiché l'integrale di dbeta è già implementato in R attraverso pbeta. Prova questo:

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