概率论：蒙特卡洛算积分（c++语言描述）

· · 学习·文化课

蒙特卡洛计算积分

原始题目

\theta = \int_{-3}^{3}\frac{1}{\sqrt{2 \pi}} e^{\frac{(x-1)^2}{2}} \theta = \int_{-\infin}^{\infin}x^{\frac{5}{3}} e^{\frac{(x-1)^2}{2}}

巨坑

参考知乎：为什么很多程序无法计算负数的立方根？

第二题如果直接求x^\frac{3}{5}行不通，至少c++里面有问题， 本代码中使用的方法是当x<0时取正数来计算，最后再减 。

代码示例

#define TEXT
#include <chrono>
#include <cmath>
#include <iostream>
#include <random>
using std::cout;
using std::endl;
using namespace std::chrono;
const unsigned int N = 1e7;
bool use1 = 1;                  // use the first function
const unsigned int seed = 520;  // random number seed

double f1(double x) {
  x--;
  double ans = x * x / 2;
  return exp(-ans);
}

double f2(double x) {
  double ans = pow((x > 0 ? x : -x), 5.0 / 3.0);
  return ans * f1(x);
}

double f(double x) {
  if (use1) {
    return f1(x);
  }
  return f2(x);
}

unsigned int work(double left, double right, double height) {
  // 初始化随机数生成器引擎
  std::mt19937 gen(seed);  // 使用Mersenne Twister算法作为随机数生成器

  // 定义一个均匀分布
  std::uniform_real_distribution<double> disx(left, right);
  std::uniform_real_distribution<double> disy(0.0, height);
  unsigned int sum = 0;
  for (int i = 0; i < N; i++) {
    double x = disx(gen);
    double y = disy(gen);
    if (f(x) > y) {
      sum++;
      if (f(x) > height) {
        std::cerr << "Warning" << endl;
      }
    }
  }
  return sum;
}

void MainCalculate1() {
  const double Con1 = sqrt(M_PI * 2.0);
  double height = 1.1, left = -3, right = 3;
  long double ans = work(left, right, height);
  ans = ans / N;
  ans = ans * (right - left) * height;
  cout << "ans=" << ans / Con1 << endl;
}

void MainCalculate2() {
  double height = 2, left = -5, right = 5;
  int ans1 = work(0, right, height);
  int ans2 = work(left, 0, height);
  long double ans = (ans1 - ans2);
  ans = ans / N * (right)*height;
  cout << "ans=" << ans << endl;
}

void text() {  // Check the parameters for soundness
  if (!use1) {
    cout << "f2(5)=" << f(5) << " f2(-5)=" << f(-5) << endl;
  }

  double ma, mi;
  for (double i = -5; i <= 5; i += 0.01) {
    ma = std::max(f(i), ma);
    mi = std::min(f(i), mi);
  }
  cout << "max=" << ma << " min=" << mi << endl;
}

int main() {
#ifdef TEXT
  text();
#endif
  auto start = system_clock::now();  // for calculate run time
  if (use1) {                        // run Main function
    MainCalculate1();
  } else {
    MainCalculate2();
  }

  auto end = system_clock::now();
  auto duration = duration_cast<microseconds>(end - start);
  cout << "used "
       << double(duration.count()) * microseconds::period::num /
              microseconds::period::den
       << " s" << endl;
  return 0;
}

运行结果

f1

max=1 min=0
ans=0.977044
used 0.929164 s

f2

f2(5)=0.00490449 f2(-5)=2.22664e-07
max=1.94433 min=0
ans=3.69945
used 3.05184 s

mathmatica验证

第一个积分

第二个积分

评价是不如用mathmatica直接秒了