ch1.R

# 環境の初期化 ----------------------------------------------------------

rm(list = ls())

## ----ch1Demo--------------------------------------------------------------------------------

## 設定と準備
rho <- 0.5 # 母相関の設定
n <- 25 # サンプルサイズ
iter <- 1000 # 反復回数

# Demo 1 ----------------------------------------------------------


# 結果を格納するオブジェクト
r <- rep(0, iter)

## シミュレーション
set.seed(123)
for (i in 1:iter) {
  Y1 <- rnorm(n, 0, 1)
  Y2 <- Y1 * rho + rnorm(n, 0, (1 - rho^2)^0.5)
  r[i] <- cor(Y1, Y2)
}

## 結果
r |> hist(xlim = c(-0.2, 1))


# Demo2 -----------------------------------------------------------

## 設定と準備
rho <- 0.5 # 母相関の設定
n <- 100 # サンプルサイズ
iter <- 1000 # 反復回数

# 結果を格納するオブジェクト
r <- rep(0, iter)

## シミュレーション
set.seed(123)
for (i in 1:iter) {
  Y1 <- rnorm(n, 0, 1)
  Y2 <- Y1 * rho + rnorm(n, 0, (1 - rho^2)^0.5)
  r[i] <- cor(Y1, Y2)
}

## 結果
r |> hist(xlim = c(-0.2, 1))


## 設定と準備
rho <- 0
n <- 2500
iter <- 10000
# 結果を格納するオブジェクト
p <- rep(0, iter)

# Demo3 -----------------------------------------------------------

## シミュレーション
set.seed(123)
for (i in 1:iter) {
  Y1 <- rnorm(n, 0, 1)
  Y2 <- Y1 * rho + rnorm(n, 0, (1 - rho^2)^0.5)
  p[i] <- cor.test(Y1, Y2)$p.value
}

## 結果
ifelse(p < 0.05, 1, 0) |> mean()


# Demo 4 ----------------------------------------------------------

## 設定と準備
rho <- 0
N <- 2500
iter <- 10000
set.seed(123)
p <- rep(0, iter)

## シミュレーション本体
for (i in 1:iter) {
  Y1 <- rnorm(n, 0, 1)
  Y2 <- Y1 * rho + rnorm(n, 0, (1 - rho^2)^0.5)
  p[i] <- cor.test(Y1, Y2)$p.value
}

## 結果
ifelse(p < 0.05, 1, 0) |> mean()


# Demo 5 ----------------------------------------------------------
## 実行に時間がかかります。ご注意ください。
## 設定と準備
rho <- 0
n <- 25
iter <- 10000
alpha <- 0.05
# 結果を格納するオブジェクト
p <- rep(0, iter)

## シミュレーション
set.seed(123)
for (i in 1:iter) {
  # 最初のデータ
  Y1 <- rnorm(n, 0, 1)
  Y2 <- Y1 * rho + rnorm(n, 0, (1 - rho^2)^0.5)
  p[i] <- cor.test(Y1, Y2)$p.value
  # データ追加
  count <- 0
  ## p値が5%を下回るか，データが当初の倍になるまで増やし続ける
  while (p[i] >= alpha && count < n * 2) {
    # 有意じゃなかったとき，それぞれの変数に1つデータを追l加
    Y1_add <- rnorm(1, 0, 1)
    Y1 <- c(Y1, Y1_add)
    Y2 <- c(Y2, Y1_add * rho + rnorm(1, 0, (1 - rho^2)^0.5))
    p[i] <- cor.test(Y1, Y2)$p.value
    count <- count + 1
  }
}

## 結果
ifelse(p < 0.05, 1, 0) |> mean()