蒙特卡洛模拟求圆周率
算法思路
代码的基本思想
是利用蒙特卡洛方法（Monte Carlo method）来估计圆周率π \piπ。蒙特卡洛方法是一种以概率统计为基础的数值计算方法，通过随机采样得到结果或近似值。在这个程序中，我们生成了一个以原点为中心、半径为r rr的圆。然后我们模拟投掷n nn个随机点，如果该点在圆内，则代表当前的点可以绕圆周作一次完整的旋转；而若该点在圆外，则表示该时刻没有绕圆周运动。最后，估计出的圆内点数p pp与总点数n nn之比再乘以4 44，即可得到圆周率的一个估计值。
代码主要包括以下几个部分：
- monte_carlo_pi函数
该函数在半径为r的单位圆内随机投掷n个点，并返回投掷到圆内的点的数量。投掷的过程实际就是根据均匀分布生成二元组的点(x,y), 然后判断该点是否在圆内，即是否满足x 2 + y 2 ⩽ r 2 x^2+y^2\leqslant r^2x2+y2⩽r2.
- 主函数main
该函数使用monte_carlo_pi函数来模拟采样，各采样用不同的随机种子。然后，将得到的圆内点数之和除以总采样数量，得到圆面积与正方形面积之比 p n \frac{p}{n}np 的近似值 p ^ \hat{p}p^. 因为圆面积为π r 2 4 \frac{\pi r^2}{4}4πr2，正方形面积为r 2 r^2r2，所以 p n ≈ π 4 \frac{p}{n} \approx \frac{\pi}{4}np≈4π.
最后再乘以4即可得到近似的π \piπ的值，并根据样本标准差和置信区间计算估计偏差。
这个代码的用处是用蒙特卡罗方法来估计圆周率。该方法可以在很短的时间内得到较为精确的结果，在数值计算中经常被使用。
完整代码
#include <stdio.h> #include <stdlib.h> #include <math.h> #include <time.h> /* 该函数模拟了在一个半径为r的圆内投掷n个点的情况，并返回样本中落在圆内的点的数量 */ int monte_carlo_pi(int n, float r) { int i; float x, y, distance; int p = 0; /* 落在圆内的点的数量 */ for (i = 0; i < n; i++) { /* 对每个点进行均匀随机分布，即随机生成一个二维点 (x, y)，且每个点被生成的概率相等 */ x = (float) rand() / RAND_MAX * 2 * r - r; /* 在[-r, r]范围内随机生成x坐标 */ y = (float) rand() / RAND_MAX * 2 * r - r; /* 在[-r, r]范围内随机生成y坐标 */ /* 判断该点是否在圆内 */ distance = sqrt(x*x + y*y); /* 计算点到原点的距离 */ if (distance <= r) { /* 如果距离小于等于r，则表示该点在圆内*/ p++; /* 投掷的点数加1 */ } } return p; } /* 主函数，对输入的参数进行处理，调用monte_carlo_pi函数模拟多次投掷情况，并输出结果 */ int main(void) { int i; /* 循环计数器 */ int n = 1000000; /* 每次模拟中投掷的点的数量 */ float r = 1.0; /* 圆的半径 */ double mu = 0, s = 0, p_mean, p_stddev, err; /* p的均值、p的方差、样本均值、样本标准差、误差程度 */ srand(time(NULL)); /* 设置随机数种子 */ /* 利用不同的随机种子进行10次模拟 */ for (i = 0; i < 10; i++) { /* 进行一次投掷，并累计落在圆内的点的数量 */ int p = monte_carlo_pi(n, r); /* 落在圆内的点的数量 */ mu += p; /* 各次投掷所得到的p的和，后面再除以10即为均值 */ s += p * p; /* 各次投掷所得到的p的平方的和，后面再除以10即为方差 */ } /* 计算得到的p的平均值和标准差 */ mu /= (n * 10); /* 根据大数定理，采样数量越多，样本均值越接近总体均值*/ // 计算样本标准差 s = s / (n * 10 - 1); // 样本方差：S^2=(∑(Xi-X)²)/(N-1) p_stddev = sqrt(s); // 样本标准差： S = √S^2 err = 1.96 * p_stddev / sqrt(n * 10); /* 根据95%置信区间计算误差程度，1.96为正态分布的值 */ /* 输出结果 */ printf("样本数量：%d\n", n * 10); /* 总共采样了多少个样本 */ printf("圆周率的估计值：%lf\n", 4 * mu); /* 各次投掷所得到的p的均值的均值 */ printf("误差程度：+-%lf\n", 4 * err); /* 置信区间为平均值加减误差程度 */ return 0; /* 程序正常结束*/ }
运行结果：
样本数量：10000000 圆周率的估计值：3.141101 误差程度：+-0.045014
