使用 gauss-hermit 求积法的积分逼近:mvQuad 包
Integral approximation using gauss-hermit quadrature method: mvQuad package
我想近似对应于 $E(XY)$ 的积分,其中 X 和 Y 是独立的,X~N(0.5,1) 和 Y~N(0.5,1),使用 gauss-hermit 正交使用mvQuad 包。
由于两个随机变量呈正态分布,近似值应该接近0.25。
但我得到了不同的结果。
> require(mvQuad)
> nw <- createNIGrid(dim=2, type="GHe", level=c(10,10))
> m=c(0.5,0.5)
> c=matrix(c(1,0,0,1),nrow = 2,byrow = F)
> rescale(nw, m = m, C = c, dec.type = 0)
>
> myFun2d <- function(x){
+ (0.5+sqrt(2)*x[,1])*(0.5+sqrt(2)*x[,2])/pi
+ }
>
> quadrature(myFun2d, grid = nw)
[1] 58.71479
>
谁能帮我弄清楚我做错了什么?
谢谢
您应该在函数 myFun2d
中为 X
和 Y
使用 PDF 的乘积
myFun2d <- function(x) {
x[, 1] * exp(-0.5 * (x[, 1]-0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2]-0.5)^2) / sqrt(2 * pi)
}
你会看到
> library(mvQuad)
> # create grid
> nw <- createNIGrid(dim = 2, type = "GHe", level = c(10, 10))
> m <- c(0.5, 0.5)
> c <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = F)
> rescale(nw, m = m, C = c, dec.type = 0)
> # define the integrand
> myFun2d <- function(x) {
+ x[, 1] * exp(-0.5 * (x[, 1] - 0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2] - 0.5)^2) / .... [TRUNCATED]
> # compute the approximated value of the integral
> (A <- quadrature(myFun2d, grid = nw))
[1] 0.25
我想近似对应于 $E(XY)$ 的积分,其中 X 和 Y 是独立的,X~N(0.5,1) 和 Y~N(0.5,1),使用 gauss-hermit 正交使用mvQuad 包。
由于两个随机变量呈正态分布,近似值应该接近0.25。
但我得到了不同的结果。
> require(mvQuad)
> nw <- createNIGrid(dim=2, type="GHe", level=c(10,10))
> m=c(0.5,0.5)
> c=matrix(c(1,0,0,1),nrow = 2,byrow = F)
> rescale(nw, m = m, C = c, dec.type = 0)
>
> myFun2d <- function(x){
+ (0.5+sqrt(2)*x[,1])*(0.5+sqrt(2)*x[,2])/pi
+ }
>
> quadrature(myFun2d, grid = nw)
[1] 58.71479
>
谁能帮我弄清楚我做错了什么?
谢谢
您应该在函数 myFun2d
X
和 Y
使用 PDF 的乘积
myFun2d <- function(x) {
x[, 1] * exp(-0.5 * (x[, 1]-0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2]-0.5)^2) / sqrt(2 * pi)
}
你会看到
> library(mvQuad)
> # create grid
> nw <- createNIGrid(dim = 2, type = "GHe", level = c(10, 10))
> m <- c(0.5, 0.5)
> c <- matrix(c(1, 0, 0, 1), nrow = 2, byrow = F)
> rescale(nw, m = m, C = c, dec.type = 0)
> # define the integrand
> myFun2d <- function(x) {
+ x[, 1] * exp(-0.5 * (x[, 1] - 0.5)^2) / sqrt(2 * pi) * x[, 2] * exp(-0.5 * (x[, 2] - 0.5)^2) / .... [TRUNCATED]
> # compute the approximated value of the integral
> (A <- quadrature(myFun2d, grid = nw))
[1] 0.25