蒙特卡罗方法计算圆周率

前几天读到了一篇网志：蒙特卡罗方法入门，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.VEVb.com/kodoyang/p/MonteCarloMethod_PI.html

雨木阳子

2015年8月9日