✨蒙特卡羅方法

蒙特·卡羅方法（Monte Carlo method），也稱統計模擬方法，是二十世紀四十年代中期由於科學技術的發展和電子計算機的發明，而被提出的一種以概率統計理論為指導的一類非常重要的數值計算方法。是指使用隨機數（或更常見的偽隨機數）來解決很多計算問題的方法。

蒙特卡羅方法是一種計算方法。原理是通過大量隨機樣本，去了解一個系統，進而得到所要計算的值。

它非常強大和靈活，又相當簡單易懂，很容易實現。對於許多問題來說，它往往是最簡單的計算方法，有時甚至是唯一可行的方法。

它誕生於上個世紀40年代美國的"曼哈頓計划"，名字來源於賭城蒙特卡羅，象征概率。

✨蒙特卡羅計算π

正方形內部有一個相切的圓，它們的面積之比是π/4。

現在，在這個正方形內部，隨機產生10000個點（即10000個坐標對 (x, y)），計算它們與中心點的距離，從而判斷是否落在圓的內部。

如果這些點均勻分布，那么圓內的點應該占到所有點的 π/4，因此將這個比值乘以4，就是π的值。通過 R語言腳本隨機模擬30000個點，π的估算值與真實值相差0.07%。

# estimating pi using Monte Carlo methods. This code won't run perfectly for you, because I 
# didn't save everything I did, like adding vectors to a dataframe named Pi. 
# Just the important stuff is here. 

center <- c(2.5, 2.5)                      # the center of the circle
radius <- 2.5
distanceFromCenter <- function(a) {
  sqrt(sum((center - a) ^ 2))
}

# let's define a 5 by 5 square matrix
points <- c(0,0, 0,5, 5,5, 5,0)
square <- matrix(points, nrow = 4, ncol = 2, byrow = TRUE)


# now all I need to do is make matrix A a matrix of real simulated points.
n <- 100                         # number of points
A <- matrix(runif(n*2, min=0, max=5), nrow = n, ncol = 2, byrow = T)   # random sampling occurs here!!!

## An alternate way to generate randoms, with faster convergence
# library(randtoolbox)
# A <- matrix(5*halton(n*2), nrow = n, ncol = 2, byrow = T)

# here's how you'll test if it's in the circle. 
b <- apply(A, 1, distanceFromCenter)           
# d <- subset(b, b < radius)                         # if you know a better way to do this part, email me.
# num <- length(d) / length(b)
num <- mean(b < radius)

piVec <- c()

for (i in 1:2000) {
  n <- i
  A <- matrix(runif(n*2, min=0, max=5), nrow = n, ncol = 2, byrow = T)
  b <- apply(A, 1, distanceFromCenter)
  d <- subset(b, b < radius)
  num <- length(d) / length(b)
  piVec[i] = num*4
}

library(data.table)
Pi <- data.frame(piVec)
Pi <- data.table(Pi)
Pi[, ind := seq(0, 1999)]
Pi[, error := abs(pi - piVec)]
Pi <- data.frame(Pi)
names(Pi) <- c("guess", "ind", "error")



##### Graphing the error

# note - if you want this part to work for you, you'll have to create the appropriate data frame from the piVec vector.
library(ggplot2)

ggplot(Pi, aes(x = ind, y = error)) +
  geom_line(colour = '#388E8E') + 
  ggtitle("Error") + 
  xlab("Sample Size") + 
  ylab("Error")

✨代碼實現

Loop實現

import random
import time

# 投擲點數總和
DARTS = 1000 * 1000
# 落在圓內的點數
hits = 0

# 開始計時
start = time.perf_counter()

for i in range(1, DARTS + 1):
    x, y = random.random(), random.random()
    dist = pow(x ** 2 + y ** 2, 0.5)
    if dist <= 1.0:
        hits = hits + 1
    else:
        pass
# 計算PI
pi = 4 * (hits / DARTS)

# 結束計時
end = time.perf_counter()

print("圓周率: {}".format(pi))
print("使用循環計算耗時: {:.5f}s".format(end - start))

Numpy性能優化

import random
import time
import numpy as np


def calPIbyLoop():
    # 投擲點數總和
    DARTS = 1000 * 1000
    # 落在圓內的點數
    hits = 0

    # 開始計時
    start = time.perf_counter()

    for i in range(1, DARTS + 1):
        x, y = random.random(), random.random()
        dist = pow(x ** 2 + y ** 2, 0.5)
        if dist <= 1.0:
            hits = hits + 1
        else:
            pass
    # 計算PI
    pi = 4 * (hits / DARTS)

    # 結束計時
    end = time.perf_counter()

    print("圓周率: {}".format(pi))
    print("使用循環計算耗時: {:.5f}s".format(end - start))


def calPIbyNumpy():
    # 投擲點數總和
    DARTS = 1000 * 1000
    # 落在圓內的點數
    hits = 0

    # 開始計時
    start = time.perf_counter()

    # 使用numpy的隨機函數生成2行DART列的矩陣: points, 代表投擲DART次點的坐標，
    # 第0行代表x軸坐標，第1行代表y軸坐標。
    points = np.random.rand(2, DARTS)
    # print(points)

    # 矩陣運算代替循環
    # numpy.where(condition[, x, y])
    # Return elements chosen from x or y depending on condition.
    hits = np.sum(np.where(((points[0] ** 2 + points[1] ** 2) ** 0.5) < 1, 1, 0))

    # 計算PI
    pi = 4 * (hits / DARTS)

    # 結束計時
    end = time.perf_counter()

    print("圓周率: {}".format(pi))
    print("使用Numpy計算耗時: {:.5f}s".format(end - start))


if __name__ == "__main__":
    calPIbyLoop()
    calPIbyNumpy()