前幾天讀到了一篇網(wǎng)志：蒙特卡羅方法入門，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);    }    }

運行結(jié)果如下：

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);    }    }

多次運行結(jié)果如下：

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;    }    }

多次運行結(jié)果如下：

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);    }}

運行結(jié)果如下：

3.148658436

參考這個級數(shù)：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);    }}

運行結(jié)果如下：

3.141592643589326

為了避免計算時間超過十秒鐘，很隨意的減小了樣本值。

【方形中的所有像素計算】中一共計算10^8次，當(dāng)在【方形中的隨機(jī)像素計算】中也計算相同的次數(shù)時，就會陷入等待。

猜測原因是獲取隨機(jī)數(shù)的時候浪費了很多時間，也可能是循環(huán)的次數(shù)太多消耗時間。

【方形中的隨機(jī)像素求平均值】中巴10^8分成了計算10^6共10^2次求平均，計算的值顯然比【方形中的所有像素計算】更接近真實的PI。

【方體中的所有小方體計算】可能是進(jìn)入了三層循環(huán)消耗了時間，當(dāng)RADIUS為1000時，實際上計算了10^9次，然而偏差很大。

【萊布尼茨公式計算】采用了公式 pi/4 =1 - 1/3 + 1/5 - 1/7 + 1/9 ...

其中，【方形中的隨機(jī)像素計算】和【方形中的隨機(jī)像素求平均值】采用了隨機(jī)的方式，屬于蒙特卡羅方法。

使用蒙特卡羅方法，可以在樣本數(shù)很小的時候也可以得到最接近真實值的值，它比全部計算更節(jié)省計算資源。

請點擊下方的『關(guān)注我』關(guān)注我！

本文地址：http://m.survivalescaperooms.com/kodoyang/p/MonteCarloMethod_PI.html

雨木陽子

2015年8月9日