definite integration with monte Carlo does not have the acceptable error

user_13519156 0
Viewed 88

I wish to calculate the integration of e^(- landa x) cosx from zero to positive infinity. The exact solution is landa/((landa^2)+1). The example I'm trying to solve forces me to have a maximum error of 0.01. Here, the approach I am taking is that I first generate a random number from uniform distribution function [0,1] and then transform it via (-1/landa ln (x)) so, each variable will now have the probability of (landa * e ^(-landa x)). The thing which I do not understand is that when I am increasing N from 1 million to be 100 million, the error changes in the following manner, which of course does not fulfill the criteria of the problem, and the weird thing is that from N of 1000000 to N of 10000000, the error is increasing. The error versus N is :

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

Here is my code:

#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; 

}
montecarlo c++
2 Answers

What @Damien wrote is good wrt the fact that you have to remove exponent part from function, you already using it for sampling.

Where he is wrong, is the fact that you have to remove whole NORMALIZED exponential PDF

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

which means that integrated function shall pick up λ in the denominator. In other words, code with cos(x) will only work for λ=1

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

In the expression above first part is PDF used for sampling, and second part is what you're trying to compute mean value.

This is fixed code, tested on 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;
}
user_40446960

You made a confusion: once you have generated x with the proper distribution, then you have to integrate cos(x) simply.

    sum = sum + cos(x);

Output:

1000000         -0.0010428
10000000        0.000105266
100000000       -2.08618e-05

Correction

As @Severin Pappadeux correctly mentioned in a comment, and in their answer, I was fast a little bit in checking a result and forgot a normalisation factor 1/landa. It has no incidence for landa = 1 of course.

There are different ways to explain this factor.

One way is the following: if you make the variable change u = exp(-landa*x), then you get

dx = -1/(landa*u) du

This landa factor in the denominator is missing in the final integration sum...

user_105533410
Related