﻿rm(list = ls())  # ワークスペースをクリア
a <- 100; b <- 30
n <- 100000  # 試行回数
set.seed(1)  # 乱数のシード値を固定
x <- runif(n, 0, a)  # U(0,a)に従う乱数をn個生成
y <- runif(n, 0, a)  # U(0,a)に従う乱数をn個生成
u <- pmin(x, y); v <- pmax(x, y)
z <- runif(n, x, a)  # U(x,a)に従う乱数をn個生成
# (1) 確率を表示 (※モンテカルロ法で求めた値，理論値の順)
sum(abs(x-y)<=b)/n ; 1 - (a-b)^2/a^2
# (2) 確率を表示 (※モンテカルロ法で求めた値，理論値の順)
sum(u<a/2 & v<u+a/2 & a/2<v)/n ; 1/4
# (3) 確率を表示 (※モンテカルロ法で求めた値，理論値の順)
sum(x>b & a-z>b)/n ; (a-2*b)/a + (b/a)*log(b/(a-b))