﻿rm(list = ls())  # ワークスペースをクリア
n <- 1000  # 試行回数
set.seed(1)  # 乱数のシード値を固定
# Pi=(xi[i], eta[i])
xi <- runif(n)  # U(0,1)に従う乱数をn個生成し，xiに代入
eta <- runif(n)  # U(0,1)に従う乱数をn個生成し，etaに代入
# 変数nに対しNn(D)を返す関数
ND <- function(n) { 
  s <- 0
  for (i in 1:n) {
    if (xi[i]^2 + eta[i]^2 <= 1) { s <- s + 1 }
  }
  return(s)
}
Z <- numeric(n)  # Nn(D)/nを格納するベクトル
for (i in 1:n) { 
  Z[i] <- ND(i) / i 
}
split.screen(1:2)
screen(1)
plot(Z, type = "l", xlab = "n", ylab = "Nn(D)/n")
# グラフの凡例を右上に追加
legend("topright", legend = c("Nn(D)/n", "mu"), lty = 1:2)
f <- function(x) { pi / 4 }  # 理論値μを返す定数関数
vf <- Vectorize(f)  # スカラー関数をベクトル化
curve(vf, add = TRUE, lty = 2)
screen(2)
# Pのプロット
plot(xi, eta, xlab = "x", ylab = "y")
circle <- function(x) { sqrt(1 - x^2) }  # x^2 + y^2 = 1 のグラフ
vcircle <- Vectorize(circle)  # スカラー関数をベクトル化
curve(vcircle, add = TRUE)