蒙特卡洛算法

蒙特卡洛算法（Monte Carlo Method）更准确的说是一种方法，是一种统计学的方法，是一种模拟的方法。

蒙特卡洛模拟是二战期间，为了解决原子弹研制工作中，裂变物资的中子随机扩散问题，美国数学家冯 诺依曼和乌拉姆等提出的一种统计方法，代号：蒙特卡洛（摩纳哥二战时期非常著名的赌城）。

• 通过大量随机样本，去了解一个系统，进而得到所要计算的值。

一、蒙特卡洛方法求π值

• 圆的面积 = PI * R * R
• 方形面积 = (2R) * (2R) = 4 * R * R
• PI = 4 * 圆 / 方

视图层

import javax.swing.*;
import java.awt.*;

public class AlgoFrame extends JFrame {

    private int canvasWidth;
    private int canvasHeight;


    public AlgoFrame(String title, int canvasWidth, int canvasHeight){
        super(title);

        this.canvasWidth = canvasWidth;
        this.canvasHeight = canvasHeight;

        Algocanvas canvas = new Algocanvas(); // 画布
        setContentPane(canvas);
        pack();

        setResizable(false);
        setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        setVisible(true);
    }

    public int getCanvasWidth() {
        return canvasWidth;
    }

    public int getCanvasHeight() {
        return canvasHeight;
    }

    // TODO: 设置自己的数据
    private MonteCarloPiData data;
    public void render(MonteCarloPiData data){
        this.data = data;
        repaint();
    }

    // 内部类自定义画布
    private class Algocanvas extends JPanel{

        public Algocanvas() {
            // 双缓存
            super(true);  // 默认为true
        }

        @Override
        public void paintComponent(Graphics g) {
            super.paintComponent(g);

            Graphics2D g2d = (Graphics2D)g;

            // 抗锯齿
            RenderingHints hints = new RenderingHints(
                    RenderingHints.KEY_ANTIALIASING,
                    RenderingHints.VALUE_ANTIALIAS_ON );
            g2d.addRenderingHints(hints);

            // 具体绘制
            // TODO: 绘制自己的数据data
           Circle circle = data.getCircle();
           AlgoVisHelper.setStrokeWidth(g2d, 3);
           AlgoVisHelper.setColor(g2d, AlgoVisHelper.Blue);
           AlgoVisHelper.strokeCircle(g2d, circle.getX(), circle.getY(), circle.getR());

            for (int i = 0; i < data.getPointsNumber(); i++) {
                Point p = data.getPoint(i);
                if(circle.contain(p))
                    AlgoVisHelper.setColor(g2d, AlgoVisHelper.Red);
                else
                    AlgoVisHelper.setColor(g2d, AlgoVisHelper.Green);

                AlgoVisHelper.fillCircle(g2d, p.x, p.y, 2);
            }
        }

        @Override
        public Dimension getPreferredSize() {
            return new Dimension(canvasWidth, canvasHeight);
        }
    }
}

控制层

import java.awt.*;

public class AlgoVisualizer {

    public static int DELAY = 100;  // 更新延迟的时间

    private MonteCarloPiData data;
    private AlgoFrame frame;
    private int N; 

    public AlgoVisualizer(int sceneWidth, int sceneHeight, int N){

        if(sceneWidth != sceneHeight)
            throw new IllegalArgumentException("This is Demo must be run in a square.");
        // 初始化数据
        this.N = N;
        Circle circle = new Circle(sceneWidth/2, sceneHeight/2, sceneWidth/2);
        data = new MonteCarloPiData(circle);

        // 初始化视图
        EventQueue.invokeLater(() -> {
            frame = new AlgoFrame("Calculate Pi", sceneWidth, sceneHeight);
            frame.setLocation(70, 70);
            new Thread( () -> {
                run();
            }).start();
        });
    }

    // 动画逻辑
    private void run(){
        for (int i = 0; i < N; i++) {

            if (i % 100 == 0) {  // 每100次进行一次输出
                frame.render(data);
                AlgoVisHelper.pasue(DELAY);
                System.out.println(data.estimatePi());
            }

            int x = (int) (Math.random() * frame.getCanvasWidth());
            int y = (int) (Math.random() * frame.getCanvasHeight());
            data.addPoint(new Point(x, y));
        }
    }

    // 测试
    public static void main(String[] args) {
        int sceneWidth = 500;
        int sceneHeight = 500;
        int N = 100000;

        AlgoVisualizer visualizer = new AlgoVisualizer(sceneWidth, sceneHeight, N);
    }
}

数据层

• Circle类

import java.awt.*;

public class Circle {

    public int x, y;
    private int r;
    public boolean isFilled = false;

    public Circle(int x, int y, int r){
        this.x = x;
        this.y = y;
        this.r = r;
    }

    public int getX(){return x;}
    public int getY(){return y;}
    public int getR(){return r;}

    public boolean contain(Point p){
        return (x - p.x) * (x - p.x) + (y - p.y) * (y - p.y) <= r * r;
    }
}

• MonteCarloPiData类

import java.awt.*;
import java.util.LinkedList;

public class MonteCarloPiData {
    private Circle circle;
    private int insideCircle = 0;
    private LinkedList<Point> points;

    public MonteCarloPiData(Circle circle){
        this.circle = circle;
        points = new LinkedList<>();
    }

    public Circle getCircle(){
        return circle;
    }

    public Point getPoint(int i){
        if (i < 0 || i >= points.size())
            throw new IllegalArgumentException("out of bound in getPoint.");
        return points.get(i);
    }

    public int getPointsNumber(){
        return points.size();
    }

    public void addPoint(Point p){
        points.add(p);
        if(circle.contain(p))
            insideCircle++;
    }

    public double estimatePi(){
        if(points.size() == 0)
            return 0.0;

        int circleArea = insideCircle;
        int squareArea = points.size();
        return (double)circleArea * 4 / squareArea;
    }
}

二、三门问题（Monty Hall Probelm）

出自美国的电视游戏节目Let’s Make a Deal。问题名字来自该节目的主持人蒙提·霍尔（Monty Hall）。参赛者会看见三扇关闭了的门，其中一扇的后面有一辆汽车，选中后面有车的那扇门可赢得该汽车，另外两扇门后面则各藏有一只山羊。当参赛者选定了一扇门，但未去开启它的时候，节目主持人开启剩下两扇门的其中一扇，露出其中一只山羊。主持人其后会问参赛者要不要换另一扇仍然关上的门。问题是：换另一扇门会否增加参赛者赢得汽车的机率？

• 换门 ==> 中奖率：2/3
• 不换门 ==> 中奖率：1/3

蒙特卡洛方法模拟

public class ThreeGatesExperiment {

    private int N;

    public ThreeGatesExperiment(int N){
        if(N <= 0)
            throw new IllegalArgumentException("N must be larger than 0.");

        this.N = N;
    }

    public void run(boolean changeDoor){
        int wins = 0;
        for (int i = 0; i < N; i++)
            if(play(changeDoor))
                wins++;
        System.out.println(changeDoor ? "change" : "Not Change");
        System.out.println("Winning rate:" + (double)wins/N);
    }

    private boolean play(boolean changeDoor){
        // Door 0, 1, 2
        int prizeDoor = (int)(Math.random() * 3);  // 玩家选择的门
        int playerChoice = (int)(Math.random() * 3);  // 奖品所在的门

        if(playerChoice == prizeDoor)
            return changeDoor ? false : true;
        else
            return changeDoor ? true : false;
    }

    public static void main(String[] args) {

        int N = 10000000;
        ThreeGatesExperiment exp = new ThreeGatesExperiment(N);
        exp.run(true);
        System.out.println();
        exp.run(false);
    }
}

计算结果：

change
Winning rate:0.6668623

Not Change
Winning rate:0.3333437

三、中奖？

在游戏里，有一种宝箱，打开这个宝箱获得传奇武器的概率是20%，现在你打开5个这样的宝箱，获得传奇武器的概率是多少？

蒙特卡洛方法模拟

public class WinningPrize {

    private double chance;
    private int playTime;
    private int N;

    public WinningPrize(double chance, int playTime, int N) {
        if(chance < 0.0 || chance > 1.0)
            throw new IllegalArgumentException("chance must be between 0 and 1.");
        if(playTime <= 0 || N <= 0)
            throw new IllegalArgumentException("palyTime or N must be larger 5.");

        this.chance = chance;
        this.playTime = playTime;
        this.N = N;
    }
    
    public void run(){
        int wins = 0;
        for (int i = 0; i < N; i++)
            if(play())
                wins++;

        System.out.println("Winning rate: " + (double)wins/N);
    }

    private boolean play(){
        for (int i = 0; i < playTime; i++)
            if(Math.random() < chance)
                return true;
         return false;
    }

    public static void main(String[] args) {

        double chance = 0.2;
        int playTime = 5;
        int N = 10000000;

        WinningPrize exp = new WinningPrize(chance, playTime, N);
        exp.run();
    }
}

计算结果：

Winning rate: 0.6721177

理论计算值：

1 - 0.8**5 = 0.67232

注： 该文为慕课网课程整理的笔记，若侵权请告知，立删。