蒙地卡羅法求 PI

December 3, 2021

蒙地卡羅為摩洛哥王國之首都,該國位於法國與義大利國境,以賭博聞名。蒙地卡羅法是一種隨機演算,是指使用隨機亂數來解決計算問題的方法。

例如,以亂數散佈點配合面積公式,可以求得近似的 PI,是簡單的蒙地卡羅法運用。隨機演算雖然會有精確度等方面的疑慮,然而有些需求可能無最佳解,或者需要耗費大量運算才能得到理想解答,若隨機演算可以快速求得某個可接受的結果,就可能採取這類演算。

在圍棋上贏得人機之戰的 AlphaGo,採用的決策過程之一,就有蒙特卡羅樹搜尋(Monte Carlo Tree Search,MCTS),而 MCTS 存在目的是為了解決搜尋空間過於巨大,難以窮舉全部子樹的問題,MCTS 子節點展開後隨機進行遊戲,根據勝負結果來更新沿路至根節點的資訊,作為下次搜尋選取節點的評估之用。

解法思路

假設圓半徑為 1,四分之一圓面積就是 PI,能包含此四分之一圓的正方形面積就為 1,如果隨意地在正方形中散佈點,這些點有些會落於四分之一圓內,假設散佈了 n 個點,在圓內的點有 c 個,依比例來算,就會得到以下的公式:

蒙地卡羅法求 PI

要判斷產生的點是否落於圓內,可以令亂數產生 X 與 Y 兩個數值,如果 X² + Y² 小於 1 就是落在圓內。

程式實作

#include <stdio.h> 
#include <stdlib.h> 
#include <time.h> 

#define N 50001

int main(void) { 
    srand(time(NULL)); 
    
    int sum = 0;
    
    int i;
    for(i = 1; i < N; i++) { 
        double x = (double) rand() / RAND_MAX; 
        double y = (double) rand() / RAND_MAX; 
        if((x * x + y * y) < 1) {
            sum++; 
        }
    } 

    printf("PI = %f\n", (double) 4 * sum / (N - 1));
    
    return 0; 
}  
import static java.lang.Math.*;
public class MonteCarlo {
    public static void main(String[] args) {
        final int N = 50001;
        int sum = 0; 
        for(int i = 1; i < N; i++) if(pow(random(), 2) + pow(random(), 2) < 1) {
            sum++;
        }
        System.out.printf("PI = %f%n", 4.0 * sum / (N - 1));  
    }
}
from random import random
N = 50001
print("PI =", 4 * len([1 for i in range(1, N) 
        if random() ** 2 + random() ** 2 < 1]) / (N - 1))
import java.lang.Math._
val N = 50000
printf("PI = %f%n", 4.0 * (for(i <- 1 to N
   if(pow(random(), 2) + pow(random(), 2) < 1)) yield 1).size / N)
N = 50000
print "PI = ", 4.0 * (1...N).map { 
    rand ** 2 + rand ** 2 < 1 ? 1 : 0 }.reduce(:+) / N
Array.prototype.reduce = function(init, f) {
    var value = init;
    for(var i = 0; i < this.length; i++) {
        value = f(value, this[i]);
    }
    return value;
};

function range(n) {
    var r = [];
    for(var i = 0; i < n; i++) {
        r[i] = i;
    }
    return r;
}

var n = 50000;
print(4 * range(n).map(function() {
    var x = Math.random();
    var y = Math.random();
    return x * x + y * y < 1 ? 1 : 0;
}).reduce(0, function(ac, elem) {
    return ac + elem;
}) / n);
import System.Random

rand gen n= take n \$ randomRs (0.0, 1.0) gen::[Float]

main = do
    gen1 <- getStdGen
    gen2 <- newStdGen
    let n = 50000
        ic = length [1 | (x, y) <- 
                 zip (rand gen1 n) (rand gen2 n), (x ** 2 + y ** 2) < 1]
    print (4 * fromIntegral ic / fromIntegral n)
from '/lib/math' import random
from '/lib/math' import pow

N = 50001
(sum = iterate(0, N)
         .select(_ -> pow(random(), 2) + pow(random(), 2) < 1)
         .length())

println('PI = {0}'.format(4 * sum / (N - 1)))