﻿rm(list = ls())  # ワークスペースをクリア
n <- 100000  # 試行回数
set.seed(1)  # 乱数のシード値を固定
X <- rmultinom(n, 10, c(1, 1, 1, 1, 1, 1))  # 多項分布に従う乱数をn個生成
X1 <- X[1, ]  # 1の目が出る回数
X2 <- X[2, ]  # 2の目が出る回数
Y <- X[1, ] + X[3, ] + X[5, ]  # 奇数の目が出る回数
Z <- X[2, ] + X[4, ] + X[6, ]  # 偶数の目が出る回数
# 期待値E(X1+X2^2)を表示 (※モンテカルロ法で求めた値，理論値の順)
mean(X1 + X2 * X2) ; 35 / 6
# 相関係数ρ(X1,X2)を表示 (※モンテカルロ法で求めた値，理論値の順)
cor(X1, X2) ; - 1 / 5
# 相関係数ρ(X1,Y)を表示 (※モンテカルロ法で求めた値，理論値の順)
cor(X1, Y) ; sqrt(5) / 5
# 相関係数ρ(X2,Y)を表示 (※モンテカルロ法で求めた値，理論値の順)
cor(X2, Y) ; -sqrt(5) / 5
# 相関係数ρ(Y,Z)を表示 (※モンテカルロ法で求めた値，理論値の順)
cor(Y, Z) ; - 1
# 確率P(X2=1,Y=3)を表示 (※モンテカルロ法で求めた値，理論値の順)
sum(X2 == 1 & Y == 3) / n ; 35 / 1458