如何在 Julia 中使用 vegas 计算具有无限边界的高维多重积分
How to compute a high dimensional multiple integral with infinite bounds using vegas in Julia
我正在尝试使用 Julia(>1400 维)计算高维积分。因此,我尝试使用函数 vegas 来执行此操作,因为它大概可以计算高维积分。然而,vegas 假设积分域是 [0,1]^n,但我的积分超过了 R^n。 vegas 的 documentation 建议更改变量,但我无法让它在多个维度上工作。
因此,如果我在 Julia 中输入以下二维积分:
using LinearAlgebra, Cuba
multinormal(x,μ,Σ) = det(2*π*Σ)^(-1/2)*exp(-1/2*(x-μ)'*inv(Σ)*(x-μ))
vegas((x,f)->f=multinormal(x,[0;0],[1 0;0 1]),2)
我得到结果
Component:
1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Number of subregions: 0
Note: The desired accuracy was reached
假设积分超过 [0,1]^2。
尝试在 [0,infinity)^2 上计算相同的积分,我尝试将建议的变量更改 here 为
vegas((x,f)->f=multinormal(x./(1 .- x),[0;0],[1 0;0 1])./(1 .-x).^2,2)
这给了我结果
Component:
1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Number of subregions: 0
Note: The desired accuracy was reached
但结果应该是 0.5 而不是 0。
如何使用 vegas 计算具有无限积分极限的多元正态分布的积分?
如果您使用 Quadrature.jl,它将自动为您执行必要的变量更改。然后你只需使用 [-Inf,Inf]
边界。查看测试中的示例:
https://github.com/SciML/Quadrature.jl/blob/master/test/inf_integral_tests.jl
我最终只是按照这里的建议对积分进行了近似:https://stats.stackexchange.com/questions/228687/approximation-expectation-integral
例如,为了计算二维多元正态分布变量 x 的期望值 E,我做了:
using Distributions, LinearAlgebra
sampleSize = 1000;
dist = MvNormal([0;0],I);
x = rand(dist,sampleSize);
E = 1/sampleSize*sum([x[:,r] for r in 1:sampleSize])
这种方法的便利之处在于它对高维度(在我的例子中 x 的维度 >1400)非常有效。
我正在尝试使用 Julia(>1400 维)计算高维积分。因此,我尝试使用函数 vegas 来执行此操作,因为它大概可以计算高维积分。然而,vegas 假设积分域是 [0,1]^n,但我的积分超过了 R^n。 vegas 的 documentation 建议更改变量,但我无法让它在多个维度上工作。
因此,如果我在 Julia 中输入以下二维积分:
using LinearAlgebra, Cuba
multinormal(x,μ,Σ) = det(2*π*Σ)^(-1/2)*exp(-1/2*(x-μ)'*inv(Σ)*(x-μ))
vegas((x,f)->f=multinormal(x,[0;0],[1 0;0 1]),2)
我得到结果
Component:
1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Number of subregions: 0
Note: The desired accuracy was reached
假设积分超过 [0,1]^2。
尝试在 [0,infinity)^2 上计算相同的积分,我尝试将建议的变量更改 here 为
vegas((x,f)->f=multinormal(x./(1 .- x),[0;0],[1 0;0 1])./(1 .-x).^2,2)
这给了我结果
Component:
1: 0.0 ± 7.025180405943273e-18 (prob.: -999.0)
Integrand evaluations: 1000
Number of subregions: 0
Note: The desired accuracy was reached
但结果应该是 0.5 而不是 0。
如何使用 vegas 计算具有无限积分极限的多元正态分布的积分?
如果您使用 Quadrature.jl,它将自动为您执行必要的变量更改。然后你只需使用 [-Inf,Inf]
边界。查看测试中的示例:
https://github.com/SciML/Quadrature.jl/blob/master/test/inf_integral_tests.jl
我最终只是按照这里的建议对积分进行了近似:https://stats.stackexchange.com/questions/228687/approximation-expectation-integral
例如,为了计算二维多元正态分布变量 x 的期望值 E,我做了:
using Distributions, LinearAlgebra
sampleSize = 1000;
dist = MvNormal([0;0],I);
x = rand(dist,sampleSize);
E = 1/sampleSize*sum([x[:,r] for r in 1:sampleSize])
这种方法的便利之处在于它对高维度(在我的例子中 x 的维度 >1400)非常有效。