モンテカルロシミュレーションの基礎
乱数によるサンプリングを繰り返して、円周率などの値を推定します。
「モンテカルロシミュレーションの基礎」はCoddyKit上の無料R Academyレッスンです。 これはレッスン3/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
Monte Carlo法とは
Monte Carloシミュレーションは、解析的に計算することが難しい量を、繰り返しランダムサンプリングして推定する方法です。モナコのカジノにちなんで名付けられ、金融、物理学、統計学、AIにおける積分、最適化、推論などに使われます。
# Core idea: approximate deterministic quantities
# using random sampling and the Law of Large Numbers
# Example: estimate probability that sum of two dice > 8
set.seed(42)
n <- 10000
die1 <- sample(1:6, n, replace = TRUE)
die2 <- sample(1:6, n, replace = TRUE)
total <- die1 + die2
# Monte Carlo estimate
mc_estimate <- mean(total > 8)
cat('MC estimate P(sum > 8):', mc_estimate, '\n')
# Exact probability
exact <- sum(outer(1:6, 1:6, '+') > 8) / 36
cat('Exact probability:', exact)ランダムな点による円周率の推定
Monte Carlo法の古典的なデモでは、単位正方形内に点をランダムに配置します。単位円の内部に入る点の割合からπ/4を推定できます。n → ∞のとき、推定値はπに収束します。
set.seed(42)
n <- 100000
# Random points in [-1, 1] x [-1, 1] unit square
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
# Check if inside unit circle: x^2 + y^2 <= 1
inside <- (x^2 + y^2) <= 1
# Pi estimate: fraction inside * area of square
pi_estimate <- 4 * mean(inside)
cat('Pi estimate:', pi_estimate, '\n')
cat('True pi:', pi, '\n')
cat('Error:', abs(pi_estimate - pi))期待値の推定
Monte Carlo積分では、X₁,...,Xₙをiidな抽出値として、E[f(X)] ≈ (1/n) Σ f(Xᵢ)とします。aからbまでのg(x)の積分では、xをUniform(a,b)からサンプリングし、E[g(X)] ≈ (b-a) * mean(g(samples))と推定します。
set.seed(42)
# Estimate integral of sin(x) from 0 to pi
# Exact value = 2
n <- 50000
x_samples <- runif(n, min = 0, max = pi)
integrand_vals <- sin(x_samples)
# MC estimate = (b-a) * mean(f(x))
mc_integral <- pi * mean(integrand_vals) # (pi - 0) * mean
cat('MC estimate of integral:', mc_integral, '\n')
cat('Exact value: 2\n')
cat('Error:', abs(mc_integral - 2), '\n')
# Standard error of the estimate
se <- pi * sd(integrand_vals) / sqrt(n)
cat('Standard error:', se)繰り返しシミュレーションのためのreplicate()
replicate(n, expr)は、シミュレーションをn回実行して結果を収集するR組み込み関数です。単純なシミュレーションではforループより簡潔に書け、ベクトルまたは行列を返します。
set.seed(42)
# Simulate sample mean of 30 N(0,1) draws
# Repeat 10000 times to study sampling distribution
sample_means <- replicate(10000, {
x <- rnorm(30) # sample of 30
mean(x) # compute mean
})
# Central Limit Theorem: sample mean ~ N(0, 1/sqrt(30))
mean(sample_means) # ~0
sd(sample_means) # ~1/sqrt(30) = 0.183
# 95% CI width
diff(quantile(sample_means, c(0.025, 0.975)))
# Compare to theoretical
2 * 1.96 / sqrt(30)大数の法則のデモ
大数の法則によれば、nが増えると標本平均は真の平均に収束します。Rでこの収束を観察すると、Monte Carlo法が機能する理由と、推定値が安定する速さを理解できます。
set.seed(42)
# Rolling a fair die: true mean = 3.5
n_max <- 10000
rolls <- sample(1:6, n_max, replace = TRUE)
cumulative_means <- cumsum(rolls) / seq_along(rolls)
# Show convergence at different sample sizes
ns <- c(10, 100, 1000, 5000, 10000)
results <- data.frame(
n = ns,
mean = cumulative_means[ns],
error = abs(cumulative_means[ns] - 3.5)
)
print(results)
# Error decreases as n increases分散減少：反対変量
一様乱数サンプルuとその補数(1-u)を反対変量のペアとして使用します。これらは負の相関を持つため、関数評価回数を増やさずに推定量の分散を最大50%減らせます。
set.seed(42)
n <- 1000
# Standard MC: estimate E[exp(U)] where U~Uniform(0,1)
# True value = e - 1 = 1.718282
u <- runif(n)
mc_std <- mean(exp(u))
# Antithetic variates: use u AND 1-u
u_anti <- runif(n/2)
mc_anti <- mean((exp(u_anti) + exp(1 - u_anti)) / 2)
cat('True value:', exp(1) - 1, '\n')
cat('Standard MC:', mc_std, '\n')
cat('Antithetic MC:', mc_anti, '\n')
# Variance comparison
var_std <- var(exp(runif(10000)))
var_anti <- var((exp(runif(5000)) + exp(1 - runif(5000)))/2)
cat('Variance ratio (anti/std):', var_anti/var_std)制御変量
制御変量とは、期待値が既知で、対象の量と相関する関数です。これを適切な倍率で引くことで分散を減らします。典型的には、平均が既知の相関したg(X)を使ってE[f(X)]を推定します。
set.seed(42)
n <- 5000
# Estimate E[exp(U)] where U~Uniform(0,1)
# True: e - 1 = 1.71828
# Control variate: g(U) = U, E[U] = 0.5
u <- runif(n)
f_vals <- exp(u) # target
g_vals <- u # control variate
# Optimal coefficient c = -Cov(f,g)/Var(g)
c_star <- -cov(f_vals, g_vals) / var(g_vals)
# Control variate estimator
mc_cv <- mean(f_vals + c_star * (g_vals - 0.5))
cat('Standard MC:', mean(f_vals), '\n')
cat('Control variate MC:', mc_cv, '\n')
cat('True value:', exp(1) - 1)
# Variance reduction factor
var(f_vals) / var(f_vals + c_star * (g_vals - 0.5))確率推定へのMonte Carlo法の利用
Monte Carlo法は、解析的に扱うことが難しい複雑な確率の推定に適しています。確率過程を何度もシミュレーションし、事象が発生した回数の割合を計算します。
set.seed(123)
# Birthday problem: P(at least 2 people share birthday)
# in a group of n people
birtday_collision <- function(n_people) {
birthdays <- sample(1:365, n_people, replace = TRUE)
length(birthdays) != length(unique(birthdays))
}
# Estimate for groups of size 10, 23, 50
sizes <- c(10, 23, 50)
for (sz in sizes) {
p <- mean(replicate(5000, birt_day_collision <- {
bd <- sample(1:365, sz, replace = TRUE)
length(bd) != length(unique(bd))
}))
cat('n =', sz, ': P(collision) ~', round(p, 3), '\n')
}幾何ブラウン運動のシミュレーション
株価は、しばしば幾何ブラウン運動としてモデル化されます。S(t+dt) = S(t) * exp((μ - σ²/2)dt + σ√dt * Z)で、Z~N(0,1)です。Monte Carlo法により、オプション価格評価のための価格経路を生成できます。
set.seed(42)
S0 <- 100 # initial price
mu <- 0.05 # annual drift
sigma <- 0.2 # annual volatility
T <- 1 # 1 year
n_steps <- 252 # daily steps
dt <- T / n_steps
# Simulate one price path
Z <- rnorm(n_steps)
log_returns <- (mu - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z
price_path <- S0 * exp(cumsum(log_returns))
cat('Final price:', round(price_path[n_steps], 2), '\n')
cat('Min price:', round(min(price_path), 2), '\n')
cat('Max price:', round(max(price_path), 2))Monte Carlo法によるオプション価格評価
Monte Carlo法でヨーロピアン・コールオプションを評価します。満期時の株価を多数シミュレーションし、ペイオフmax(S_T - K, 0)を計算してから、平均ペイオフをe^(-rT)で割り引きます。
set.seed(42)
S0 <- 100; K <- 105; r <- 0.05; sigma <- 0.2; T <- 1
n_sim <- 50000
# Final stock prices under risk-neutral measure
Z <- rnorm(n_sim)
ST <- S0 * exp((r - 0.5 * sigma^2) * T + sigma * sqrt(T) * Z)
# Call option payoff
payoff <- pmax(ST - K, 0)
# Discounted expected payoff
call_price <- exp(-r * T) * mean(payoff)
se <- exp(-r * T) * sd(payoff) / sqrt(n_sim)
cat('Call price:', round(call_price, 4), '\n')
cat('95% CI: [', round(call_price - 1.96*se, 4),
',', round(call_price + 1.96*se, 4), ']')Monte Carlo法の収束速度
Monte Carlo法はO(1/√n)の速度で収束します。誤差を半分にするには、サンプル数を4倍にする必要があります。Monte Carlo推定値の標準誤差はσ/√nです。ここでσは被積分関数の標準偏差です。
set.seed(42)
# Demonstrate MC convergence for pi estimation
estimate_pi <- function(n) {
x <- runif(n, -1, 1)
y <- runif(n, -1, 1)
4 * mean(x^2 + y^2 <= 1)
}
# Sample sizes: powers of 10
ns <- 10^(1:5)
estimates <- sapply(ns, function(n) {
set.seed(42)
estimate_pi(n)
})
data.frame(
n = ns,
pi_estimate = round(estimates, 5),
error = round(abs(estimates - pi), 5)
)クイックチェック
Monte Carloシミュレーションの基本について、理解度を確認しましょう。
まとめ：Monte Carlo法の基本
重要なポイント：Monte Carlo法では、多数のランダムサンプルの平均を使って量を推定します。繰り返しシミュレーションにはreplicate()を使用します。大数の法則により収束が保証されます。誤差は1/√nに比例するため、誤差を半分にするにはサンプル数を4倍にします。分散減少手法（反対変量、制御変量）を使うと、サンプル数を増やさずに効率を高められます。応用例には積分、確率推定、オプション価格評価、シミュレーションがあります。
set.seed(42)
# Monte Carlo template:
# 1. Define simulation function
simulate_once <- function() {
x <- runif(1, -1, 1)
y <- runif(1, -1, 1)
(x^2 + y^2) <= 1
}
# 2. Replicate many times
n <- 10000
results <- replicate(n, simulate_once())
# 3. Estimate quantity of interest
pi_mc <- 4 * mean(results)
# 4. Quantify uncertainty
se <- 4 * sd(results) / sqrt(n)
cat('Pi:', pi_mc, '+/-', round(1.96*se, 4))よくある質問
「モンテカルロシミュレーションの基礎」レッスンは無料ですか？
はい。「モンテカルロシミュレーションの基礎」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「モンテカルロシミュレーションの基礎」で何を学びますか？
乱数によるサンプリングを繰り返して、円周率などの値を推定します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか？
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン3/4です。
「モンテカルロシミュレーションの基礎」レッスンにはどのくらい時間がかかりますか？
ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか？
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。
このコースのすべてのレッスン
- set.seed() と再現性
- 確率分布からの乱数生成
- モンテカルロシミュレーションの基礎
- R でのブートストラップ再標本化