蒙特卡罗方法计算圆周率
前几天读到了一篇网志：蒙特卡罗方法入门，http://www.ruanyifeng.com/blog/2015/07/monte-carlo-method.html
其中介绍了用概率计算圆周率的方法，所以就用程序做了以下尝试。
作为常量的PI值的近似在Math.PI中为3.141592653589793。
Ⅰ.方形中的所有像素计算
package yumu.probability.montecarlo;public class CalculatePI { private static final int RADIUS = 10000; public static void main(String[] args) { int circle = 0; for(int i = 0; i < RADIUS; ++i){ for(int j = 0; j < RADIUS; ++j){ if(i*i + j*j < RADIUS*RADIUS){ ++circle; } } } double quarterPI = (double)circle / (RADIUS * RADIUS); System.out.println(quarterPI * 4); } }
运行结果如下：
3.14199016
Ⅱ.方形中的随机像素计算
package yumu.probability.montecarlo;public class RandomPI { private static final int QUANTITY = 10000000; public static void main(String[] args) { double x, y; int circle = 0; for(long i = 0; i < QUANTITY; ++i){ x = Math.random(); y = Math.random(); if(x*x + y*y < 1){ ++circle; } } double quarterPI = (double)circle / QUANTITY; System.out.println(quarterPI * 4); } }
多次运行结果如下：
3.14113443.1422143.1410763.1407648
Ⅲ.方形中的随机像素求平均值
package yumu.probability.montecarlo;public class RandomTwicePI { private static final int COUNT = 100; private static final int QUANTITY = 1000000; public static void main(String[] args) { double sum = 0; for(int i = 0; i < COUNT; ++ i){ sum += randomPI(); } double pi = sum / COUNT; System.out.println(pi); } private static double randomPI(){ double x, y; int circle = 0; for(long i = 0; i < QUANTITY; ++i){ x = Math.random(); y = Math.random(); if(x*x + y*y < 1){ ++circle; } } double quarterPI = (double)circle / QUANTITY; return quarterPI*4; } }
多次运行结果如下：
3.14175815999999933.1414950800000013.1415940399999998
Ⅳ.方体中的所有小方体计算
package yumu.probability.montecarlo;public class CalculateCubePI { private static final int RADIUS = 1000; public static void main(String[] args) { int sphere = 0; for(int i = 0; i < RADIUS; ++i){ for(int j = 0; j < RADIUS; ++j){ for(int k = 0; k < RADIUS; ++k){ if(i*i + j*j + k*k < RADIUS*RADIUS){ ++sphere; } } } } double oneInSixOfPI = (double)sphere / (RADIUS * RADIUS * RADIUS); System.out.println(oneInSixOfPI * 6); }}
运行结果如下：
3.148658436
Ⅴ.莱布尼茨公式计算
参考这个级数：Leibniz formula for π, https://en.wikipedia.org/wiki/Leibniz_formula_for_%CF%80
package yumu.probability.montecarlo;public class LeibnizSeriesPI { private static final int COUNT = 100000000; public static void main(String[] args) { double quarterPI = 1; for(int i = 1; i < COUNT; ++i){ double temp = i%2==0 ? 1 : -1; temp /= i*2 + 1; quarterPI += temp; } System.out.println(quarterPI * 4); }}
运行结果如下：
3.141592643589326
Ⅵ.总结
为了避免计算时间超过十秒钟，很随意的减小了样本值。
【方形中的所有像素计算】中一共计算10^8次，当在【方形中的随机像素计算】中也计算相同的次数时，就会陷入等待。
猜测原因是获取随机数的时候浪费了很多时间，也可能是循环的次数太多消耗时间。
【方形中的随机像素求平均值】中巴10^8分成了计算10^6共10^2次求平均，计算的值显然比【方形中的所有像素计算】更接近真实的PI。
【方体中的所有小方体计算】可能是进入了三层循环消耗了时间，当RADIUS为1000时，实际上计算了10^9次，然而偏差很大。
【莱布尼茨公式计算】采用了公式 pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 ...
其中，【方形中的随机像素计算】和【方形中的随机像素求平均值】采用了随机的方式，属于蒙特卡罗方法。
使用蒙特卡罗方法，可以在样本数很小的时候也可以得到最接近真实值的值，它比全部计算更节省计算资源。
请点击下方的『关注我』关注我！
本文地址：http://www.cnblogs.com/kodoyang/p/MonteCarloMethod_PI.html
雨木阳子
2015年8月9日