﻿rm(list = ls())  # ワークスペースをクリア
f <- function(x){ sin(x) }  # 被積分関数f(x)
# モンテカルロ法による積分計算
monte_carlo <- function(n, a) {
  sigma <- 0
  set.seed(3)  # 乱数のシード値を固定
  for (i in 1:n) {
    Z <- runif(1, 0, a)  # U(0,a)に従う乱数を1個生成し，Zに代入
    sigma <- sigma + f(Z)
  }
  return(sigma / n)
}
# プロット関数
mc_plot <- function(n, a) {
  y <- numeric(n)
  for(i in 1:n) {
    # ベクトルyにmonte_carlo(i, a) (i=1,2,...,n)を順に代入
    y[i] <- monte_carlo(i, a) 
  }
  # (i,y[i])をプロットし，直線で結ぶ
  plot(y, main = paste("n_max =", n), 
       xlab = "n", ylab = "Monte Carlo integration", 
       type = "l") 
}
a <- pi  # 円周率π
mc_plot(1000, a)
integral_val <- 2  # ∫[0,a]f(x)dx=∫[0,π]sin(x)dxの理論値
# MC積分の系列が理論値に収束する様子を描画
abline(h = integral_val / a, col = "blue")  # h=線の縦軸位置