exp(-x^2/2) 从 x=-infinity 到 x=+infinity 的 Monte Carlo 积分

Question

我要整合

f(x) = exp(-x^2/2)

从 x=-infinity 到 x=+infinity

使用蒙特卡洛方法。我使用函数 randn() 为函数 f(x_i) = exp(-x_i^2/2) 生成所有 x_i 我想积分以计算 f([x_1,..x_n]) 的平均值。我的问题是，结果取决于我为边界 x1 和 x2 选择的值（见下文）。通过增加 x1 和 x2 的值，我的结果与实际值相差甚远。其实通过增加x1和x2结果应该会越来越好。

有人看到我的错误吗？

这是我的 Matlab 代码

clear all;
b=10;                 % border
x1 = -b;              % left border
x2 = b;               % right border
n = 10^6;             % number of random numbers
x = randn(n,1);
f = ones(n,1);
g = exp(-(x.^2)/2); 
F = ((x2-x1)/n)*f'*g;

正确的值应该是 ~2.5066。

谢谢

Answer 1

试试这个：

clear all;
b=10;                 % border
x1 = -b;              % left border
x2 = b;               % right border
n = 10^6;             % number of random numbers

x = sort(abs(x1 - x2) * rand(n,1) + x1);
f = exp(-x.^2/2);
F = trapz(x,f)

F =

    2.5066

Answer 2

好的，让我们开始编写 MC 集成的一般案例：

I = S f(x) * p(x) dx, x in [a...b]

S 这里是整数符号。

通常p(x)是归一化概率密度函数，f(x)要积分，算法很简单：

• 将累加器s 设置为零
• 开始循环N 事件
• 从p(x)随机抽样x
• 给定x，计算f(x) 并添加到累加器中
• 如果没有完成则返回开始循环
• 如果完成，将累加器除以N 并返回

在最简单的教科书案例中

I = S f(x) dx, x in [a...b]

这意味着 PDF 等于均匀分布的一个

p(x) = 1/(b-a)

但你要求和的实际上是(b-a)*f(x)，因为你的积分现在看起来像

I = S (b-a)*f(x) 1/(b-a) dx, x in [a...b]

一般来说，如果f(x) 和p(x) 都可以用作PDF，那么您可以选择将f(x) 集成到p(x) 上，还是将p(x) 集成到f(x) 上。没有不同！ （嗯，除了计算时间）

所以，回到特定的积分（我相信它等于\sqrt{2\pi}）

I = S exp(-x^2/2) dx, x in [-infinity...infinity]

您可以使用更传统的方法，例如 @Agriculturist 并编写它

I = S exp(-x^2/2)*(2a) 1/(2a) dx, x in [-a...a]

并在 [-a...a] 区间内从 U(0,1) 中采样 x，并为每个 x 计算 exp() 并对其进行平均并得到结果

据我了解，您想将exp() 用作PDF，因此您的积分看起来像

I = S D * exp(-x^2/2)/D dx, x in [-infinity...infinity]

要归一化的PDF，因此它应包含归一化因子D，它完全等于高斯积分的\sqrt{2 \pi}。

现在f(x) 只是一个等于D 的常数。它不依赖于x。这意味着您应该为每个采样的x 添加一个 CONSTANT 值D 到累加器。运行N 样本后， 在累加器中，您将完全拥有N*D。要找到平均值，您将除以N，结果您将得到完美的D，即\sqrt{2 \pi}，反过来，它是 2.5066.

写任何matlab太生疏了，无论如何新年快乐