definite integration with monte Carlo does not have the acceptable error

Question

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; 

}

Answer 1

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