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

Question

我要整合

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

从 x=-infinity 到 x=+infinity

通过使用 Monte Carlo 方法。我使用函数 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 和 return 它

在最简单的教科书案例中你有

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 添加到累加器的常数值 D。在运行 N 个样本之后，在累加器中，您将恰好 N*D。要找到均值，您将除以 N，结果您将得到完美的 D，即 \sqrt{2 \pi}，而这又是 2.5066.

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

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

Monte Carlo integration of exp(-x^2/2) from x=-infinity to x=+infinity

matlab

integration

normal-distribution

montecarlo