﻿rm(list = ls())  # ワークスペースをクリア
n <- 100000  # 試行回数
b1_denominator <- 0 ; b1_numerator <- 0  # (1)の計算で利用
b2_denominator <- 0 ; b2_numerator <- 0  # (2)の計算で利用
b3_denominator <- 0 ; b3_numerator <- 0  # (3)の計算で利用
set.seed(1)  # 乱数のシード値を固定
for(i in 1:n) {
  x <- sample(c(0,1), 1)  # 上の子のデータ (0は男, 1は女)
  y <- sample(c(0,1), 1)  # 下の子のデータ (0は男, 1は女)
  if(x == 0) {  # (1) 上の子が男
    b1_denominator <- b1_denominator + 1
    if(y == 0) {  # 下の子も男
      b1_numerator <- b1_numerator + 1
    }
  }
  if(x * y == 0) {  # (2) 少なくとも1人が男
    b2_denominator <- b2_denominator + 1
    if(x + y == 0) {  # どちらも男
      b2_numerator <- b2_numerator + 1
    }
  }
  z <- sample(c(0,1), 1)  # (3) 0なら上の子, 1なら下の子
  if((z == 0 && x == 0) || (z == 1 && y == 0)) {  # 無作為に選んだ1人が男
    b3_denominator <- b3_denominator + 1
    if(z == 0 && y == 0) {  # もう1人の子も男
      b3_numerator <- b3_numerator + 1
    }
    if(z == 1 && x == 0) {  # もう1人の子も男
      b3_numerator <- b3_numerator + 1
    }
  }
}
# (1) モンテカルロ法で求めた確率を表示 (※参考：理論値は1/2)
b1_numerator / b1_denominator
# (2) モンテカルロ法で求めた確率を表示 (※参考：理論値は1/3)
b2_numerator / b2_denominator
# (3) モンテカルロ法で求めた確率を表示 (※参考：理論値は1/2)
b3_numerator / b3_denominator