标签：

蒙特卡罗方法计算圆周率

前几天读到了一篇网志：蒙特卡罗方法入门，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.1411344
3.142214
3.141076
3.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.1417581599999993
3.141495080000001
3.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日

蒙特卡罗方法计算圆周率

标签：

原文地址：http://my.oschina.net/kodoyang/blog/494402