我尝试执行函数的蒙特卡罗积分,但采样不均匀。我需要这种方法既适用于对数尺度,也适用于极坐标的积分,因为我会将两者结合起来,并在半径中使用极坐标与对数尺度采样。
我写了一个测试脚本
- 在极坐标中积分2D高斯函数(应等于pi)
- 从10**-2到107以对数尺度积分y(x)= x(其应等于~0.5*1014)
为了测试的目的,我补充了一个统一的笛卡尔坐标为基础的Monte Carlo的工作计算。正是样本的不均匀性改变了我的结果。
import numpy as np
def function_to_integrate(x, y):
return np.exp(-x**2 - y**2)
def polar_MC(polar):
size = 100000
integral = 0.
integration_radius = 4.
if polar:
for _ in range(size):
r = np.random.random()*integration_radius
phi = np.random.random()*2.*np.pi
x = r*np.cos(phi)
y = r*np.sin(phi)
jacobian_MC_polar = 1.
integral += function_to_integrate(x, y) * jacobian_MC_polar
integral = integral * np.pi * integration_radius**2 / size
else:
for _ in range(size):
length = 2. * integration_radius
x = np.random.random()*length - length/2.
y = np.random.random()*length - length/2.
integral += function_to_integrate(x, y)
integral = integral * length**2 / size
print('POLAR: True integral should be pi ', '; MC:', integral, polar)
def log_MC(log):
size = 10000
integral = 0.
if log:
for _ in range(size):
x = np.random.uniform(-2, 7.)
jacobian_MC_log = 1.
integral += 10**x * jacobian_MC_log
else:
for _ in range(size):
x = np.random.uniform(10**-2, 10**7)
integral += x
integral = integral*10**7 / size
print('LOG: True integral should be 0.5*10**7*10**7 = 5*10**13; MC:', integral/10**13, '* 10**13', log)
polar_MC(polar=True)
polar_MC(polar=False)
log_MC(log=True)
log_MC(log=False)
字符串
我无法从极坐标和对数尺度蒙特卡罗得到正确的结果,我应该如何设置jacobian_MC才能使其工作?还是我做错了什么?
我试过使用标准的雅可比矩阵(r表示极坐标,r* np.log(10)表示对数),但这并没有帮助。
当雅可比数设为1时
POLAR: True integral should be pi ; MC: 11.041032315593327 True
POLAR: True integral should be pi ; MC: 3.108344559871783 False
LOG: True integral should be 0.5*10**7*10**7 = 5*10**13; MC: 0.48366198481209793 * 10**13 True
LOG: True integral should be 0.5*10**7*10**7 = 5*10**13; MC: 5.003437412553992 * 10**13 False
型
增加采样没有帮助,结果接近收敛。
我应该用什么概率分布来划分采样点?
1条答案
按热度按时间mxg2im7a1#
对于极积分,你把雅可比和归一化部分都搞错了
这是正确的代码,Python 3.10,Win x64
字符串
它总是打印出3.14左右的值:
型
更新
这不是边界或区域的事情。这是如何制作概率密度函数(PDF)的。
所以f(r)
f(r)= r e-r2
和积分
I = S02 pi d phi S0R dr f(r)
你想用蒙特卡洛。这意味着对φ和r进行采样。#21518;样的东西(任何东西)。)你必须有PDF(正的,正确地标准化为1等)。
让我们从对φ的积分开始
Iphi = S02pi d phi。
我将用2 π乘除它。
Iphi = 2pi S02pi d phi/(2pi)。
在集成标志下制作正确的PDF
PDF(phi)= d phi/(2 pi)
和适当的取样
phi = 2 pi random()
但是剩下的是前面的2 π,你必须把它推到前面。这就是它的工作原理-你在积分下归一化为1 PDF,无论那里有什么常数,你都要在整个采样例程之后考虑它。它不是面积也不是周长,而是PDF构造和标准化
同样适用于径向积分
IR = S0R dr = R S dr/R。
PDF(r)= dr/R,所以你可以采样r = R random(),但是你必须把R提前带到最后的计算步骤。