蒙地卡羅方法(MonteCarlomethod)，是一種以機率與統計理論出發的數值計算方法，這門課程將透過幾個範例介紹如何在R環境中用蒙地卡羅方法解決問題。
蒙地卡羅方法主要能處理兩類問題：(1).所求解問題可以轉化為某種隨機分佈的特徵數，譬如隨機事件出現的概率，或者隨機變數的期望值。(2).所求解的問題本身具有內在的隨機性，藉助電腦直接模擬這種隨機的過程，譬如分析中子在核反應堆中的傳輸過程。
範例1.蒙地卡羅積分。給定母體平均數mean=0,標準差sd=10的常態分佈，計算函數定義域介於(-10,10)之間的面積。積分結果為0.682689，可以參考這個網址http://goo.gl/T62O4W
請同學依照註解上的指示修改參數xx,yy以計算蒙地卡羅積分，完成後，請存檔並回到console輸入submit()
set.seed(1) # 給定亂數生成種子
runs <- 1000 # 模擬次數
# 請同學修改 xx 成適當的參數，以生成平均數=0, 標準差=10的常態分佈隨機變數
sims <- rnorm(runs, mean=0, sd=10)
# 請同學修改 yy 成適當的參數，以計算介於 (-10, 10)的機率總和
mc.integral <- sum(sims >= -10 & sims <= 10)/runs
stopifnot(class(mc.integral) == "numeric")
# 完成後，請存檔並回到console輸入`submit()`同學利用蒙地卡羅法所計算出來的積分結果已經被存到mc.integral這個變數，請同學輸入該變數以觀察成果
mc.integral透過wolframalpha得到積分結果為0.682689，與我們僅模擬1000次的結果仍有差距，積分結果可以參考這個網址http://goo.gl/T62O4W
試著增加模擬次數，我們透過函數與向量運算計算蒙地卡羅積分在模擬次數一千次到一百萬次的結果，完成後，請存檔並回到console輸入submit()
set.seed(1) # 給定亂數生成種子
runs <- 10^(3:6) # 模擬次數
# 這是一個給定參數 runs=模擬次數，即執行MC積分演算法的函數
mcFun01 <- function(runs){
# 生成平均數=0, 標準差=10的常態分佈隨機變數
sims <- rnorm(runs, mean=0, sd=10)
# 計算介於 (-10, 10)的機率總和
mc.integral <- sum(sims >= -10 & sims <= 10)/runs
mc.integral
}
# 請同學修改 xx 與 yy 成適當的參數，計算不同模擬次數的MC積分結果
mc.integral.seq <- sapply(X=runs, FUN=mcFun01)
stopifnot(class(mc.integral.seq) == "numeric")
stopifnot(length(mc.integral.seq) == "4")
# 完成後，請存檔並回到console輸入`submit()`請同學對runs與mc.integral.seq作圖，觀察模擬次數與積分結果的關係，請輸入以下程式碼plot(log(runs,10),mc.integral.seq,type="l")
plot(log(runs,10), mc.integral.seq, type="l", ylim=c(0.67,0.69))請同學在圖片加上真實的積分結果，請輸入以下程式碼abline(h=0.682689)
abline(h=0.682689)請問同學，增加模擬次數是不是會讓蒙地卡羅積分更準確？ Yes
範例2.計算A/B測試的p值。考慮A、B兩種網頁排版方案，比較會員註冊的成效。在同一時間內隨機分派訪客進入A、B兩種排版的網頁，其中A版本在2,981個瀏覽數中有343個註冊數；B版本在2,945個瀏覽數中有403個註冊數。A版本的轉換率為11.5%；B版本為13.7%，
練習用統計檢定的方式驗證A/B版本是否ㄧ具有統計上的顯著差異。實務上，我們會假設轉換率滿足分佈beta(k+1,n-k+1)，其中k=註冊數,n=瀏覽數，欲了解分佈假設的理論可參考這個網址http://goo.gl/yK6h3P
請同學依註解的指示生成beta分佈之隨機變數，並計算版本A大於版本B的相對比率。該比率即蒙地卡羅的近似p值。完成後，請存檔並回到console輸入submit()
set.seed(1) # 給定亂數生成種子
runs <- 100000 # 模擬次數
n_a <- 2981 # 請填入版本A的瀏覽次數
k_a <- 343 # 請填入版本A的註冊次數
n_b <- 2945 # 請填入版本B的瀏覽次數
k_b <- 403 # 請填入版本B的註冊次數
# 生成版本A與B的隨機變數
a.samples <- rbeta(runs, k_a+1, n_a-k_a+1)
b.samples <- rbeta(runs, k_b+1, n_b-k_b+1)
mc.pvalue <- sum(a.samples > b.samples)/runs
stopifnot(class(mc.pvalue) == "numeric")
# 完成後，請存檔並回到console輸入`submit()`
同學利用蒙地卡羅法所計算出來的積分結果已經被存到mc.pvalue這個變數，請同學輸入該變數以觀察成果
mc.pvalue更直覺比較AB版本轉換率差異的方式是做圖，請輸入以下程式碼hist(b.samples-a.samples)，畫出透過蒙地卡羅法生成AB版本差異隨機變數的直方圖
hist(b.samples-a.samples)再加上分隔線，請輸入以下程式碼abline(v=0,col=2)
abline(v=0, col=2)從圖中可以發現，分隔線(a-b=0)左邊的面積幾乎可略，其面積正好是mc.pvalue，P值可以解釋成AB版本沒有差異的假設下，落在拒絕域的最小機率。欲了解P值的定義請參考這個網址https://en.wikipedia.org/wiki/P-value
範例3.預測股價。我們將利用蒙地卡羅方法進行一個預測蘋果股價的實作。
在介紹這個範例前需要先了解布朗運動(Brownianmotion)，又稱Wiener過程。對於理論部分想要深入了解的同學請參考黃文璋教授的專文http://goo.gl/0yPrt6。
請依照指示生成一個布朗運動的隨機過程，step1.生成1000筆標準常態分配的隨機變數：x<-rnorm(1000,0,1)，其中這1000筆隨機變數是彼此獨立的(i.i.d.)。
x <- rnorm(1000, 0, 1)step2.建構布朗運動的隨機過程：bm<-cumsum(x)，其中cumsum是向量的累加函數。
bm <- cumsum(x)step3.畫出布朗運動的結果：plot(bm,type="l")，其中cumsum是向量的累加函數。
plot(bm, type="l")bm變數是透過彼此獨立的常態分佈隨機變數累加而成，所以它有一個明顯的特徵，就是不重複間隔的兩兩差值是獨立的，即d[s]=bm[s+1]-bm[s]與d[t]=bm[t+1]-bm[t]，在s<t時彼此獨立。
step4.利用自相關函數(Autocorrelation,一個訊號在兩次觀察之間的相關程度)，可以檢驗這個特徵，請同學輸入：acf(diff(bm))，其中diff是差值函數(後減前)，acf是自相關函數。
acf(diff(bm))從自相關的圖形中可以發現，在diff(bm)只跟自己有關(Lag=0,Cor=1)，跟其他時間無關(Lag!=0,Cor=0)。關於自相關的說明請參考https://en.wikipedia.org/wiki/Autocorrelation
step5.檢驗diff(bm)呈現標準常態分佈，請同學輸入hist(diff(bm))觀察直方圖的趨勢。
hist(diff(bm))接下來我們正是進入蘋果股價的實作，我們預先準備好一份取自Yahoofinance的蘋果歷史股價數據，請同學輸入head(AAPL)印出前6筆數據，資料包含6個欄位：每天開盤價(AAPL.Open)、最高價(AAPL.High)、最低價(AAPL.Low)、收盤價(AAPL.Close)、成交量(AAPL.Volume)、(還原權值)AAPL.Adjusted
head(AAPL)如果像自行抓取股價歷史資料的同學，可以使用quantmod這個R套件，我們提供的蘋果股價數據是用下面這段指令取得的：quantmod::getSymbols("AAPL")，關於quantmod的運用，可以參考吳牧恩教授的文章http://www.bituzi.com/2014/12/Rbacktest6mins.html
step1.請同學輸入以下指令：AAPL.price<-AAPL[,6]，以擷取還原權值AAPL.Adjusted
AAPL.price <- AAPL[,6]step2.請同學輸入以下指令：plot(AAPL.price,type="l")，畫出蘋果的歷史股價
plot(AAPL.price, type="l")從圖中可以發現，股價並非一個平穩的狀態(stationary)，改從回報率著手，我們使用logpricereturns公式，對這種回報率指標有興趣的同學請參考http://goo.gl/tB7vnJ
step3.請同學輸入以下指令：AAPL.returns<-diff(log(AAPL.price))，計算回報率
AAPL.returns <- diff(log(AAPL.price))step4.檢驗回報率的時間序列，請同學輸入以下指令：plot(AAPL.returns,type="l")
plot(AAPL.returns, type="l")回報率的時間序列有呈現水平震盪的趨勢。
step5.檢驗回報率的分佈，請同學輸入hist(AAPL.returns,breaks=40)觀察直方圖的趨勢，其中參數breaks表示切割直方圖區間的份數。
hist(AAPL.returns, breaks=40)觀察直方圖可以發現，回報率AAPL.returns呈現對稱分佈，扣除極端值後還滿像常態分佈的。
step6.檢驗回報率的自相關，請同學輸入acf(AAPL.returns)。
acf(AAPL.returns)請問同學，自相關的結果是不是顯示不同時段的回報率彼此獨立？ Yes
我們的目的是利用蒙地卡羅方法模擬蘋果的股價，進而做到股價預測。我們也發現股價與布朗運動的行為很類似。
第t期的回報率R_t=log(P_t/P_t−1)，透過轉換可得第t期的股價P_t=P_t−1*exp(R_t)。
step8.請同學依註解的指示生成一次股價模擬數據。完成後，請存檔並回到console輸入submit()
set.seed(1) # 給定亂數生成種子
N <- 30 # 模擬時間長度
mR <- mean(AAPL.returns) # 修改xxx, 計算AAPL.returns的平均數
sdR <- sd(AAPL.returns) # 修改xxx, 計算AAPL.returns的標準差
# 請執行股價的蒙地卡羅模擬
# Hint: P_t = P_t−1 * exp(R_t)
# P_t = P_0 * prod(exp(R_t))
# P_t = P_0 * exp(cumsum(R_t))
P0 <- tail(AAPL.price,1) # 設定模擬起始時間，為最後一筆股價
Rt <- rnorm(n=N, mean=mR, sd=sdR) # 修改xxx, 生成AAPL.returns的常態分佈隨機變數
Pt <- P0*exp(cumsum(Rt)) # 請同學依照Hint生成股價模擬數據
plot(Pt, type="l") # 畫出預測結果
# 請注意保持起始亂數種子 set.seed(1)，已通過syopifnot的驗證
stopifnot(length(Rt) == N)
stopifnot(class(Pt) == "numeric")
stopifnot(round(Pt[1],4) == 104.3704)
stopifnot(round(Pt[N],4) == 114.6805)
# 完成後，請存檔並回到console輸入`submit()`最後，我們展示生成多次股價模擬數據以進行預測的方法。請同學直接執行rstatistics-02-03-02.R的範例碼，看完結果後回到console輸入submit()已完成本次課程。
set.seed(1) # 給定亂數生成種子
runs <- 10 # 模擬次數
N <- 30 # 模擬時間長度
mR <- mean(AAPL.returns) # 修改xxx, 計算AAPL.returns的平均數
sdR <- sd(AAPL.returns) # 修改xxx, 計算AAPL.returns的標準差
# 請執行股價的蒙地卡羅模擬
# Hint: P_t = P_t−1 * exp(R_t)
# P_t = P_0 * prod(exp(R_t))
# P_t = P_0 * exp(cumsum(R_t))
P0 <- tail(AAPL.price,1) # 設定模擬起始時間，為最後一筆股價
# 生成 runs 筆蒙地卡羅模擬結果
tmp <- rnorm(n=N*runs, mean=mR, sd=sdR)
Rt.arr <- matrix(tmp, ncol=runs)
Pt.arr <- apply(Rt.arr, 2, function(Rt) P0*exp(cumsum(Rt)))
matplot(Pt.arr, type="l", col="grey", lty=1) # 畫出 runs 次模擬
lines(rowMeans(Pt.arr), type="l", lwd=1.5) # 畫出平均模擬結果
stopifnot(round(rowMeans(Pt.arr)[1],4) == 106.1569)
stopifnot(round(rowMeans(Pt.arr)[N],4) == 111.7722)
# 本檔案為範例頁面，請直接存檔並回到console輸入`submit()`