如何解决自定义 Python Monte Carlo 积分函数低估了多维积分
我需要创建一个自定义的 Monte Carlo 积分函数来适应使用 NumPy 的自定义多维分布对象。我需要它在每个维度上整合相同的值。它对单个维度工作正常,但在多个维度上低估,随着维度的增加,情况变得更糟。我使用这个 paper(方程 5)作为指导。我的体积 * 平均密度方程不正确吗?我的采样方法不正确吗?我真的不知道错误是什么。
import numpy as np
from scipy.stats import multivariate_normal
# Set up distribution parameters.
dim = 3
loc = np.repeat(0.,repeats=dim)
scale = np.repeat(1.,repeats=dim)
# Initialize a multivariate normal distribution.
mvn = multivariate_normal(mean=loc,cov=scale)
def mc_integrator(distribution,dim,support,size=1000,seed=0):
"""
Parameters
----------
distribution : function
A probability density function.
dim : int
The number of dimensions of the distribution.
support : list
List of the low and high values of the hypercube to integrate over.
size : int,optional
Number of samples used for estimation. The default is 1000.
seed : int,optional
A random seed for reproducibility. The default is 0.
Returns
-------
float
The estimate of the integral over the hypercube.
"""
# Set the random seed.
np.random.seed(seed)
# Separate the elements of the support.
a,b = support[0],support[1]
# Calculate the volume of the hypercube.
volume = (b-a)**dim
# Generate random samples of the appropriate shape.
samples = np.random.uniform(low=a,high=b,size=(size,dim))
# Return the estimate of the integral.
return volume*np.mean(distribution(samples))
# Set the number of samples to use for estimation.
size = 10000
# Set the low and high value over each dimension of the hypercube.
support = [-2,2]
# Print the estimate of the integral.
print(mc_integrator(mvn.pdf,size=size))
# Print the exact value of the integral.
print(mvn.cdf(np.repeat(support[1],dim))-mvn.cdf(np.repeat(support[0],dim)))
Output: 0.8523870204938726
0.9332787601629401
解决方法
John,整体看起来不错,但在我看来,您计算出的预期结果有误。我认为预期的结果应该是 (F(2) - F(-2)^3
,其中 F
是均值 0 和方差 1 的高斯 cdf。对于 F(2) - F(-2)
,我得到 erf(sqrt(2))
,大约为 0.9545,然后(F(2) - F(-2))^3
为 0.8696,与您的结果非常吻合。
我不知道 mvn.cdf
应该返回什么,但是“cdf”的概念在不止一个维度上有点可疑,所以也许你可以避开它。
关于一般的多维积分,您提到使用 Halton 序列。我认为这也是一个有趣的想法。我在计算积分方面的经验是在 1 维或 2 维中使用正交规则,在 3 到几个(5?7?我不知道)中使用低差异序列,而 MC 则不止于此。哦,还有我的建议是,在诉诸数值近似之前,要努力获得精确的结果。
我很想知道你在做什么。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。