与 monte Carlo 的定积分没有可接受的错误

Question

我想计算 e^(- landa x) cosx 从零到正无穷大的积分。精确解是landa/((landa^2)+1)。我试图解决的示例迫使我的最大误差为 0.01。在这里，我采取的方法是，我首先从均匀分布函数 [0,1] 中生成一个随机数，然后通过 (-1/landa ln (x)) 对其进行变换，因此每个变量现在的概率为 (landa * e ^(-landa x))。我不明白的是，当我将 N 从 100 万增加到 1 亿时，错误会按以下方式变化，这当然不符合问题的标准，奇怪的是从 N 1000000 到 10000000 的 N 个，误差在增加。与 N 的误差是：

N=1000000    0.0997496
N=100000000  0.0999462
N=100000000  0.0999341

这是我的代码：

#include <iostream>
#include <random>
#include <fstream>
#include <iomanip>
using namespace std;
double landa = 1;
double function(double x) {
    return (exp(-landa * x) * cos(x));
}

int main()
{
    unsigned seed = 0;
    srand(seed);
    double exact_solution = landa / (pow(landa, 2) + 1);
    const int N = 100000000;
    default_random_engine g(seed);
    uniform_real_distribution<double> distribution(0.0f, nextafter(1.0f, DBL_MAX));
    double sum = 0.0;
    double app;
    double error;
    for (int i = 0; i < N; i++) {
    double x = distribution(g);
    // transform xs 
    x = (-1.0 / landa) * log(x);
    sum = sum + function(x);
    }

    app = sum / static_cast<double> (N);
    error = exact_solution - app;
    cout << N << "\t" << error << endl; 

}

Answer 1

你弄糊涂了：一旦你用正确的分布生成了 x，那么你必须简单地集成 cos(x)。

    sum = sum + cos(x);

输出：

1000000         -0.0010428
10000000        0.000105266
100000000       -2.08618e-05

更正

正如@Severin Pappadeux 在评论中正确提到的那样，在他们的回答中，我检查结果的速度有点快，忘记了归一化因子 1/landa。 landa = 1当然没有发生。

有不同的方式来解释这个因素。

一种方法如下：如果你让变量改变u = exp(-landa*x)，那么你会得到

dx = -1/(landa*u) du

分母中的这个 landa 因子在最终积分和中缺失...

Answer 2

@Damien 写的很好，因为您必须从函数中删除指数部分，您已经在使用它进行采样。

他错的地方在于，您必须删除整个 NORMALIZED 指数 PDF

PDF(x|λ) = λ*exp(-λ x),

表示积分函数取分母中的λ。换句话说，带有 cos(x) 的代码仅适用于 λ=1

f(x) = λ*exp(-λ x) * cos(x)/λ

在上面的表达式中，第一部分是用于采样的 PDF，第二部分是您要计算的平均值。

这是固定代码，已在 Win10 x64、LLVM 11.0.1 上测试

#include <cmath>
#include <random>
#include <iostream>

using namespace std;

double lambda = 1.5;

inline double function(double x) {
    return cos(x) / lambda;
}

int main() {

    double exact_solution = lambda / (lambda*lambda + 1.0);

    const int N = 100000000;

    std::mt19937_64 g(12345671LL);

    uniform_real_distribution<double> distribution{}; // a bit raster with range [0...1)

    double sum = 0.0;
    for (int i = 0; i != N; ++i) {
        double x = distribution(g);
        x = (-1.0 / lambda) * log(1.0 - x); // to avoid domain error
        sum += function(x);
    }

    double app   = sum / static_cast<double> (N);
    double error = exact_solution - app;
    cout << N << "\t" << error << endl;
    return 0;
}

与 monte Carlo 的定积分没有可接受的错误

definite integration with monte Carlo does not have the acceptable error

c++

montecarlo