在R錯誤中的雙重集成

Question

有人可以指出我在這段代碼中出錯的地方嗎？ 我試圖解決這個積分在R：

m1 <- 70 
m2 <- 70 
s1 <- 4 
s2 <- 4.5 
rho <- 0.25 

InnerFunc <- Vectorize(function(x,y) { 

    q <- (1/(1-rho^2)) * (((x-m1)/s1)^2 -2 * rho * ((x-m1)/s1) * ((y-m2)/s2) + ((y-m2)/s2)^2) 

    Vectorize(function(y){ 
    0.009136982 * exp(-q/2) 
    }) 
    }) 

InnerIntegral <- Vectorize(function(x) { 
    Vectorize(function(y) {integrate(InnerFunc, 75, 100)$value}) 
    }) 

integrate(InnerIntegral, 75, 100) 

這導致了以下錯誤消息：

Error: Error in integrate(InnerIntegral, 75, 100) : 
evaluation of function gave a result of wrong type 

Answer 1

第一步驟是定義被積作爲的函數兩個變量如下：

integrand <- function(x, y) { 

    m1 <- 70 
    m2 <- 70 
    s1 <- 4 
    s2 <- 4.5 
    rho <- 0.25 

    q <- (1/(1-rho^2)) * (((x-m1)/s1)^2 -2 * rho * ((x-m1)/s1) * ((y-m2)/s2) + ((y-m2)/s2)^2) 
    f <- 1/(2 * pi * 4 * 4.5 * sqrt(1 - 0.0625)) 

    return (f * exp(-q/2)) 
} 

請注意，這裏不需要使用Vectorize()，因爲算術和exp()已經被矢量化了。下一步驟是定義另一個功能，其中執行上的積分y：

integrand2 <- Vectorize(function(x) { 
    ix <- function(y) integrand(x, y) 
    return (integrate(ix, 75, 100)$value) 
}) 

這裏，Vectorize()是有用的，因爲它代表了功能需要x是一個標量（或者說，向量長度1）。最後，可以執行x的集成：

integrate(integrand2, 75, 100)