C++ 蒙特卡洛求积分

1. 最简单的蒙特卡洛法求指数函数在[0,1]上的积分，面积法

const int numberOfMCSimulations = 1000000;
double lowerBound = 0.0;
double upperBound = 1.0;
// We need to draw a rectangle which, for positive functions, has a range on the y--axis from 0.0 till (at least) the maximum of the function within the rectanagle
// Since we integrate the exponential function here, which is monotonic, the maximum is simply the value of the function at the right side of the rectangle
double maxYRectangle = exp(upperBound);
double minYRectangle = 0.0;
// This variable counts the simuated MC points that end up under the curve
double pointsUnderCurve = 0.0;

// The random numbers are generated in the same way as for the Shape project
std::random_device random_dev;
std::mt19937 random_gen(random_dev());
// We sample from uniform distributions on a rectangle in the XY-plane. We have put explicitly the bounds of the rectangle
// Note that you have a uniform distribution in any finite interval, you can simply performa linear transformation to map it to any other interval (this is what we did in the first shape project)
std::uniform_real_distribution<double> distX(lowerBound, upperBound);
std::uniform_real_distribution<double> distY(minYRectangle, maxYRectangle);

for (int i = 0; i < numberOfMCSimulations; ++i)
	{
		double randX = distX(random_gen);
		double randY = distY(random_gen);
		if (exp(randX) >= randY)
		{
			++pointsUnderCurve;
		}
	}

double integralThorughMC = (upperBound - lowerBound) * (maxYRectangle - minYRectangle) * pointsUnderCurve / numberOfMCSimulations;

std::cout << "Integral_0^10 exp analytically: " << exp(upperBound) - exp(lowerBound) << std::endl;
std::cout << "Integral_0^10 exp calculated with basic MC integrator: " << integralThorughMC << std::endl;

但这个方法有个缺陷，我必须知道此函数的最值才能确定长方形的范围，所以采用改进的方法。

2. 这里假设 a=0,b=1
   

void integration_standard()
{

	const int numberOfMCSimulations = 1000000;

	double lowerBound = 0.0;
	double upperBound = 1.0;

	std::random_device random_dev;
	std::mt19937 random_gen(random_dev());
	// We sample from uniform distributions on a rectangle in the XY-plane
	std::uniform_real_distribution<double> distX(lowerBound, upperBound);

	double functionValues = 0.0;
	for (int i = 0; i < numberOfMCSimulations; ++i)
	{
		// We simply sum all function values of the uniformly distributed numbers together
		functionValues += exp(distX(random_gen));

	}


	// The approximation of the integral is simply the area of the rectangle times the average function value under the uniform distribution in the integration interval
	double integralThorughMC = (upperBound - lowerBound) * functionValues / numberOfMCSimulations;

	// For the exp function we know the analytic result
	std::cout << "Integral_0^10 exp analytically: " << exp(upperBound) - exp(lowerBound) << std::endl;
	// Now we output the MC result
	std::cout << "Integral_0^10 exp calculated with basic MC integrator: " << integralThorughMC << std::endl;
}