﻿rm(list = ls())  # ワークスペースをクリア
a <- 0.5
# 大数の強法則：m1
expXvars <- function(n) {
  set.seed(1)  # 乱数のシード値を固定
  X <- runif(n, -a, a)  # U(-a,a)に従う乱数をn個生成し，Xに代入
  return(sum(exp(X)) / n)
}
v_expXvars <- Vectorize(expXvars)  # スカラー関数をベクトル化
curve(v_expXvars, from = 0, to = 200000, 
      xlab = "n", ylab = "expXvars(n)", lty = 1)
m1 <- function(n) { (exp(a) - exp(-a)) / (2 * a) }
v_m1 <- Vectorize(m1)  # スカラー関数をベクトル化
curve(v_m1, lty = 2, add = TRUE)
legend("topright", legend = c("expXvars", "m1"), lty = 1:2)
# 大数の強法則：m2
XSQvars <- function(n) {
  set.seed(1)  # 乱数のシード値を固定
  X <- runif(n, -a, a)  # U(-a,a)に従う乱数をn個生成し，Xに代入
  return(sum(X * X) / n)
}
v_XSQvars <- Vectorize(XSQvars)  # スカラー関数をベクトル化
curve(v_XSQvars, from = 0, to = 200000, 
      xlab = "n", ylab = "XSQvars(n)", lty = 1)
m2 <- function(n) { a^2 / 3 }
v_m2 <- Vectorize(m2)  # スカラー関数をベクトル化
curve(v_m2, lty = 2, add = TRUE)
legend("topright", legend = c("XSQvars", "m2"), lty = 1:2)
# 大数の強法則：m3
XSQvars_over_expXvars <- function(n) {
  set.seed(1)  # 乱数のシード値を固定
  X <- runif(n, -a, a)  # U(-a,a)に従う乱数をn個生成し，Xに代入
  return(sum(X * X) / sum(exp(X)))
}
# スカラー関数をベクトル化
v_XSQvars_over_expXvars <- Vectorize(XSQvars_over_expXvars)
curve(v_XSQvars_over_expXvars, from = 0, to = 200000, 
      xlab = "n", ylab = "XSQvars_over_expXvars(n)", lty = 1)
m3 <- function(n) { 2 * a^3 / (3 * (exp(a) - exp(-a))) }
v_m3 <- Vectorize(m3)  # スカラー関数をベクトル化
curve(v_m3, lty = 2, add = TRUE)
legend("topright", 
       legend = c("XSQvars_over_expXvars", "m3"), lty = 1:2)