使用Monte Carlo与scipy.integrate.nquad的不同积分结果

您重新分配分布

sample = kernel.resample(size=50000)

然后计算每个采样点的概率小于边界处的概率

insample = kernel(sample) < iso

这是不正确的。考虑边界（0,100），并假设您的数据具有u =（0,0）和cov = [[100,0]，[0,100]]。点（0,50）和（50,0）在此内核中具有相同的概率，但是其中只有一个在边界内。由于两者都通过了测试，因此您采样过度。

您应该测试每个点是否sample在边界内，然后计算概率。就像是

def int1(kernel, x1, y1):
    # Sample KDE distribution                                                                                                              
    sample = kernel.resample(size=100)

    include = (sample < np.repeat([[x1],[y1]],sample.shape[1],axis=1)).all(axis=0)
    integral = include.sum() / float(sample.shape[1])
    return integral

我使用以下代码对此进行了测试

def measure(n):

    m1 = np.random.normal(size=n)
    m2 = np.random.normal(size=n)
    return m1,m2

a = scipy.stats.gaussian_kde( np.vstack(measure(1000)) )
print(int1(a,-10,-10))
print(int2(a,-10,-10))
print(int1(a,0,0))
print(int2(a,-0,-0))

产量

0.0
(4.304674927251112e-232, 4.6980863813551415e-230)
0.26
(0.25897626178338407, 1.4536217446381293e-08)

蒙特卡洛集成应该像这样工作

• 在x / y可能值的某些子集上抽样N个随机值（一致地，不是来自于您的分布）（在下面，我将其与均值相差10个SD）。
• 对于每个随机值计算内核（rand_x，rand_y）
• 计算总和并乘以（体积）/ N_samples

在代码中：

def mc_wo_sample(kernel,x1,y1,lboundx,lboundy):
    nsamples = 50000
    volume = (x1-lboundx)*(y1-lboundy)
    # generate uniform points in range                                                                                                     
    xrand = np.random.rand(nsamples,1)*(x1-lboundx) + lboundx
    yrand = np.random.rand(nsamples,1)*(y1-lboundy) + lboundy
    randvals = np.hstack((xrand,yrand)).transpose()
    print randvals.shape
    return (volume*kernel(randvals).sum())/nsamples

运行以下

   print(int1(a,-9,-9))
   print(int2(a,-9,-9))
   print(mc_wo_sample(a,-9,-9,-10,-10))
   print(int1(a,0,0))
   print(int2(a,-0,-0))
   print(mc_wo_sample(a,0,0,-10,-10))

产量

0.0
(4.012958496109042e-70, 6.7211236076277e-71)
4.08538890986e-70
0.36
(0.37101621760650216, 1.4670898180664756e-08)
0.361614657674