C++ Monte Carlo 集成:如何在不求和结果的情况下多次运行代码?
Posted
技术标签:
【中文标题】C++ Monte Carlo 集成:如何在不求和结果的情况下多次运行代码?【英文标题】:C++ Monte Carlo Integration: How to run code multiple times without summing the results? 【发布时间】:2019-05-09 07:37:06 【问题描述】:这是我在将其更改为包含多次运行之前的代码:
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <random>
#include <iomanip>
#include <vector>
using namespace std;
int main ()
double vol;
double hit;
int samples;
int i, j;
double sum;
double pt;
double actual_vol;
const double PI = 2.0*atan2(1.0,0.0);
double abs_err;
double rel_err;
random_device dev;
default_random_engine e dev() ;
uniform_real_distribution<double> u0.0,1.0;
samples = 1000000 * dim;
actual_vol = pow(PI, double(dim/2.0)) / exp(lgamma(double(dim/2.0)+1.0));
for (i = 0; i < samples; i++)
sum = 0;
for (j = 0; j < dim; j++)
pt = 2*(u(e)-0.5);
sum += pt*pt;
if (sqrt(sum) < 1)
hit += 1;
vol = ( pow(2,dim) * hit ) / samples;
abs_err = fabs( actual_vol - vol);
rel_err = abs_err / actual_vol;
cout << "Average volume of your sphere: " << setprecision(7) << vol << endl;
cout << "Actual volume: " << setprecision(7) << actual_vol << endl;
cout << "Absolute Error: " << setprecision(7) << abs_err << endl;
cout << "Relative Error: " << setprecision(7) << rel_err << endl;
我会得到如下所示的正确输出:
Average volume of your sphere: 3.140924
Actual volume: 3.141593
Absolute Error: 0.0006686536
Relative Error: 0.000212839
现在,当我更改它以便我可以调用该函数并多次运行它时,使用以下代码:
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <math.h>
#include <random>
#include <iomanip>
#include <vector>
using namespace std;
double monte_carlo (int dim)
double vol;
double hit;
int samples;
int i, j;
double sum;
double pt;
double actual_vol;
const double PI = 2.0*atan2(1.0,0.0);
double abs_err;
double rel_err;
random_device dev;
default_random_engine e dev() ;
uniform_real_distribution<double> u0.0,1.0;
samples = 1000000 * dim;
actual_vol = pow(PI, double(dim/2.0)) / exp(lgamma(double(dim/2.0)+1.0));
for (i = 0; i < samples; i++)
sum = 0;
for (j = 0; j < dim; j++)
pt = 2*(u(e)-0.5);
sum += pt*pt;
if (sqrt(sum) < 1)
hit += 1;
vol = ( pow(2,dim) * hit ) / samples;
abs_err = fabs( actual_vol - vol);
rel_err = abs_err / actual_vol;
cout << "Average volume of your sphere: " << setprecision(7) << vol << endl;
cout << "Actual volume: " << setprecision(7) << actual_vol << endl;
cout << "Absolute Error: " << setprecision(7) << abs_err << endl;
cout << "Relative Error: " << setprecision(7) << rel_err << endl;
int main (int argc, char* argv[])
int dim = 0;
int runs = 0;
int i;
dim = atoi(argv[1]);
runs = atoi(argv[2]);
for (i = 0; i < runs; i++)
monte_carlo(dim);
return 0;
我得到了这些结果,现在将以前的值与当前值相加,这不是我想要的:
Average volume of your sphere: 3.141764
Actual volume: 3.141593
Absolute Error: 0.0001713464
Relative Error: 5.454126e-05
Average volume of your sphere: 6.283674
Actual volume: 3.141593
Absolute Error: 3.142081
Relative Error: 1.000156
Average volume of your sphere: 9.427502
Actual volume: 3.141593
Absolute Error: 6.285909
Relative Error: 2.000867
Average volume of your sphere: 12.56937
Actual volume: 3.141593
Absolute Error: 9.427775
Relative Error: 3.000954
Average volume of your sphere: 15.71272
Actual volume: 3.141593
Absolute Error: 12.57113
Relative Error: 4.001515
Average volume of your sphere: 18.85378
Actual volume: 3.141593
Absolute Error: 15.71219
Relative Error: 5.001345
Average volume of your sphere: 21.99504
Actual volume: 3.141593
Absolute Error: 18.85345
Relative Error: 6.001239
您会注意到球体平均体积的第一个值约为 3.14,然后是第二个实例,现在是 6.28(或第一个实例的两倍),第三个实例是 9.42(大约是第一个实例的三倍)一)等。
它应该做的是每次运行都运行一个新的计算,并且每个计算的值都应该在 3.14 左右徘徊。如何让它停止对上一次运行的值求和?
谢谢!!
【问题讨论】:
【参考方案1】:这可能是因为您从未重新初始化变量。
您还具有强烈的“旧”C 偏见(C 标头,atoi
、fabs
... 的用法),在需要时声明变量,而且您的路径总是相似的,因为您使用具有相同种子的相同随机数生成器(默认构造)。
不过,对于您的问题:
double hit = 0;
double samples = 0;
等等。
同样对于 PI,如果您有 boost,请使用它的常数,而不是以低于它所能提供的精度重新计算它。
【讨论】:
天哪,我不敢相信我错过了那个,谢谢老兄,就是这样。 另外,强烈的 C 偏向会导致什么后果?有什么负面的吗? 由于相同的随机数生成器,如何解决路径相似的问题? 我的建议是将此函数放在一个将数字生成器作为成员变量的类中。 @KeithJF1 你的 C 偏见的后果之一(它只在 old C 中需要)是在函数顶部定义一堆变量会大大增加更容易使用未定义的值。如果你等到有一个值来初始化它之后才定义一个变量,你不能。以上是关于C++ Monte Carlo 集成:如何在不求和结果的情况下多次运行代码?的主要内容,如果未能解决你的问题,请参考以下文章
python inefficient_monte_carlo.py
Markov Chain Monte Carlo 和 Gibbs Sampling算法