自定义PythonMonte Carlo积分函数低估多维积分

Question

我需要创建自定义 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, dim, support, 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

Answer 1

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 应该是什么 return，但是“cdf”的概念在不止一个维度上有点可疑，所以也许你可以避开那个.

关于一般的多维积分，您提到使用 Halton 序列。我认为这也是一个有趣的想法。我在计算积分方面的经验是在 1 维或 2 维中使用正交规则，在 3 到几个（5？7？我不知道）中使用低差异序列，而 MC 不止于此。哦，我的建议是在求助于数值近似值之前努力获得准确的结果。

我很想知道你在做什么。

自定义PythonMonte Carlo积分函数低估多维积分

Custom Python Monte Carlo integration function is underestimating multi-dimensional integrals

numpy

probability

numerical-methods

montecarlo

numerical-integration