﻿rm(list = ls())  # ワークスペースをクリア
# Ynの漸近分布である正規分布N(μ,σ^2)のパラメータμ
mu <- function(p) { 100 * p }
# Ynの漸近分布である正規分布N(μ,σ^2)のパラメータσ
sigma <- function(n, p) { 100 * sqrt(p * (1 - p) / n) }
# 受験者数m，設問数n，各設問の得点率pを引数として
# Zn，Ynのm人分の結果をベクトルZ，Yに格納し，リスト(Z,Y)を作成
genlist_ZY <- function(m, n, p) {  
  Z <- numeric(m); Y <- numeric(m)
  # 100 * p，100 * sqrt(p * (1 - p) / n)
  mu_p = mu(p) ; sigma_n_p = sigma(n, p)
  # np，sqrt(np(1-p)) 
  np = n * p ; sqrt_npq = sqrt(n * p * (1 - p))
  set.seed(1)  # 乱数のシード値を固定
  for(i in 1:m) {
    X <- rbinom(n, 1, p) # Be(p)に従うn個の乱数X1, X2, ..., Xnを生成
    Z[i] <- (sum(X) - np) / sqrt_npq
    Y[i] <- mu_p + sigma_n_p * Z[i]
  }
  list_ZY <- list(Z, Y)  # Z,Yをリスト化
  return(list_ZY)  # リスト(Z,Y)を返す
}
# 受験者数m，設問数n，各設問の得点率pを指定
m <- 100000 ; n <- 500 ; p <- 0.6
# Zn，Ynのm人分の結果をベクトルZ，Yに格納し，リスト(Z,Y)を作成
listZY <- genlist_ZY(m, n, p)
split.screen(1:2)
# Znのヒストグラム
screen(1)
hist(listZY[[1]], main = "Histogram of Zn (n=500)", xlab = 'Zn', freq = FALSE)
curve(dnorm, add = TRUE)  # N(0,1)の密度関数を重ねて表示
# Ynのヒストグラム
screen(2)
hist(listZY[[2]], main = "Histogram of Yn (n=500)", xlab = 'Yn', freq = FALSE)
# N(μ,σ^2)の密度関数を重ねて表示
curve(dnorm(x, mu(p), sigma(n, p)), add = TRUE)  