Python实现D维蒙特卡洛积分结果与精确值不符，求错误排查

我正尝试通过蒙特卡洛积分求解如下积分：

核心思路是生成N个采样点，按如下公式计算曲线下面积：

为此我编写了如下Python代码：

import numpy as np
from sympy import symbols, integrate     

def f(x,D):                             
  return D*(x**2)

for i in range(1, 9):                    
  x = symbols('x')                      
  print("D为{}时积分的精确数学值为：{}\n".format(i, integrate(f(x,i),(x, 0,1)).evalf(2))) 

print("*************************************************************************\n")

N = 10**4
for j in range(1,9):
  ans = 0
  n_tot = N
  n_below_curve = 0
  for i in range(N):
    x0=np.random.uniform(0,1)
    y0=np.random.uniform(0,1)
    if (f(x0,j) <= y0):
      n_below_curve += 1
  ans = ( n_below_curve / n_tot ) * (1*1)
  print("D为{}时的积分结果为：{}\n".format(j, ans))

运行输出如下：

D为1时积分的精确数学值为：0.33 

D为2时积分的精确数学值为：0.67 

D为3时积分的精确数学值为：1.0 

D为4时积分的精确数学值为：1.3 

D为5时积分的精确数学值为：1.7 

D为6时积分的精确数学值为：2.0 

D为7时积分的精确数学值为：2.3 

D为8时积分的精确数学值为：2.7 

************************************************************************* 

D为1时的积分结果为：0.6635 

D为2时的积分结果为：0.4681 

D为3时的积分结果为：0.3823 

D为4时的积分结果为：0.3321 

D为5时的积分结果为：0.2978 

D为6时的积分结果为：0.269 

D为7时的积分结果为：0.252 

D为8时的积分结果为：0.2372 

对比精确值和蒙特卡洛输出结果，积分完全失效，请问代码错误出在哪里？

代码存在两个核心逻辑错误：

1. 投点判断条件写反
   你要统计的是落在曲线下方的点，也就是采样的y值小于等于函数值，即y0 <= f(x0,j)，但你写的是f(x0,j) <= y0，刚好统计了曲线上方的点，这也是D=1时结果约为0.66，刚好是1减去精确值0.33的原因。

2. y轴采样范围不匹配函数最大值
   你的被积函数f(x,D) = Dx²在积分区间[0,1]的最大值为D（当x=1时取得），但你将y的采样范围限制在了[0,1]，当D>1时，大量函数值超过了1，这部分区间的点你根本无法正确统计，自然结果完全偏离。同时计算面积时，采样矩形的面积应该是x区间长度乘以y区间长度，也就是1D，而非你写的11。

import numpy as np
from sympy import symbols, integrate     

def f(x,D):                             
  return D*(x**2)

for i in range(1, 9):                    
  x = symbols('x')                      
  print("D为{}时积分的精确数学值为：{}\n".format(i, integrate(f(x,i),(x, 0,1)).evalf(2))) 

print("*************************************************************************\n")

N = 10**4
for j in range(1,9):
  n_tot = N
  n_below_curve = 0
  y_max = j # 函数最大值为D=j
  for i in range(N):
    x0=np.random.uniform(0,1)
    y0=np.random.uniform(0,y_max) # 调整y的采样范围到0到函数最大值
    if y0 <= f(x0,j): # 修正判断条件
      n_below_curve += 1
  ans = (n_below_curve / n_tot) * (1 * y_max) # 修正采样矩形面积计算
  print("D为{}时的积分结果为：{:.2f}\n".format(j, ans))

运行修正后的代码，输出结果会和精确值基本一致。

内容的提问来源于stack exchange，提问作者J.Snowden