﻿import numpy as np  # モジュールNumpyをnpという名前でインポート
n = 100000  # 試行回数
np.random.seed(1)  # 乱数のシード値を固定
X = np.random.multinomial(10, [1/6] * 6, size = n)  # 多項分布に従う乱数
X1 = X[:, 0]  # 1の目が出る回数
X2 = X[:, 1]  # 2の目が出る回数
Y = X[:, 0] + X[:, 2] + X[:, 4]  # 奇数の目が出る回数
Z = X[:, 1] + X[:, 3] + X[:, 5]  # 偶数の目が出る回数
# 期待値E(X1+X2^2)を表示
print(np.mean(X1 + X2 * X2))  # モンテカルロ法
print(35 / 6)  # 参考：理論値
# 相関係数ρ(X1,X2)を表示
print(np.corrcoef(X1,X2)[0,1])  # モンテカルロ法
print(-1/5)  # 参考：理論値
# 相関係数ρ(X1,Y)を表示
print(np.corrcoef(X1,Y)[0,1])  # モンテカルロ法
print(np.sqrt(5) / 5)  # 参考：理論値
# 相関係数ρ(X2,Y)を表示
print(np.corrcoef(X2,Y)[0,1])  # モンテカルロ法
print(-np.sqrt(5) / 5)  # 参考：理論値
# 相関係数ρ(Y,Z)を表示
print(np.corrcoef(Y,Z)[0,1])  # モンテカルロ法
print(-1)  # 参考：理論値
# 確率P(X2=1,Y=3)を表示
print(np.sum((X2==1) & (Y==3)) / n)  # モンテカルロ法
print(35 / 1458)  # 参考：理論値