読者です 読者をやめる 読者になる 読者になる

廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

モンテカルロ法の収束のオーダーについて少し調べた

R 数値計算

動機

シミュレーションとか, とりあえず10000回くらい回しとけばいいのかな, くらいの認識しかないので, ちゃんと知っときたい.

定義;オーダー(ランダウの記号)

ふたつの関数 f: \mathbb{N} \to \mathbb{R}g: \mathbb{N} \to \mathbb{R} があるとする.

 f(N)=\mathcal{O}(g(N)) と表す(これを f\mathcal{O}(g) のオーダーであるという)とき,その意味は N \to \infty とした場合, 定数  N_0 \in \mathbb{N} および c>0 が存在して, 任意の N \ge N_0 に対して,

\displaystyle |f(N)| \le c|g(N)|

となることである.

モンテカルロ法による期待値の推定

 X_1, \ldots, X_N が独立同分布から得られたランダム標本とし, f(X) を実数上の関数とすると, 期待値  E(f(X)) は,

\displaystyle Z_{N}^{MC} = \frac{1}{N} \sum_{j=1}^{N} f(X_j)

で得られる.

Z_{N}^{MC} の平均二乗誤差は 1/N のオーダー(\mathcal{O}(1/N)), 平均二乗誤差(MSE)は 1/N のオーダー(\mathcal{O}(1/N)), 二乗平均平方根誤差(RMSE)は 1/\sqrt{N} のオーダー(\mathcal{O}(1/\sqrt{N}))である.

例;R による数値実験

100個の標準正規乱数, 1000個の標準正規乱数, 1000個の標準正規乱数から平均を求めることをそれぞれ10000回繰り返し, 精度を見る.

m_100 <- sapply(1:10000, function(i)mean(rnorm(100))) #100個
m_1000 <- sapply(1:10000, function(i)mean(rnorm(1000)))
m_10000 <- sapply(1:10000, function(i)mean(rnorm(10000)))

MSE_100 <- (m_100-0)^2
MSE_1000 <- (m_1000-0)^2
MSE_10000 <- (m_10000-0)^2

#nihongo()
boxplot(cbind("100個" = MSE_100,"1000個" =MSE_1000,"10000個"= MSE_10000), main="MSE")
abline(h=1/100,lty=2, col="red3", lwd=2)
abline(h=1/1000,lty=2, col="blue4", lwd=2)
abline(h=1/10000,lty=2, col="orange3", lwd=2)
legend("topright",legend=expression(1/100,1/1000,1/10000),
       lty=2, col=c("red3","blue4", "orange3"),lwd="2")

RMSE_100 <- sqrt((m_100-0)^2)
RMSE_1000 <- sqrt((m_1000-0)^2)
RMSE_10000 <- sqrt((m_10000-0)^2)

boxplot(cbind("100個" = RMSE_100,"1000個" =RMSE_1000,"10000個"= RMSE_10000),main="RMSE")
abline(h=sqrt(1/100),lty=2, col="red3", lwd=2)
abline(h=sqrt(1/1000),lty=2, col="blue4", lwd=2)
abline(h=sqrt(1/10000),lty=2, col="orange3", lwd=2)
legend("topright",legend=expression(sqrt(1/100), sqrt(1/1000), sqrt(1/10000)),
       lty=2, col=c("red3","blue4", "orange3"),lwd="2")

f:id:abrahamcow:20150116191416p:plain

f:id:abrahamcow:20150116192026p:plain

参考文献

ランダウの記号 - Wikipedia

(主に pp. 74-80 を参照)