In [1]:
# 確率と統計
In [1]:
# Import basic modules
import numpy as np
import scipy as sp
import pandas as pd
from pandas import Series, DataFrame
# import modules for visualization
%matplotlib inline
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns # Samuel Norman Seaborn
# show
%precision 3
# To fix random seed
np.random.seed(0)
In [6]:
# サイコロの目 = 全ての可能な根元事象を集めた標本空間S
s = [1, 2, 3, 4, 5, 6]
dice_value = np.array(s)
# サイコロの目リストから1回抽出する = 試行
np.random.choice(dice_value) # 5 ( 試行結果 = 根元事象,基本事象と呼ぶ )
# 複数回の試行も引数を与えることで可能
np.random.choice(dice_value, 2) # array([1, 4])
np.random.choice(dice_value, 3) # array([4, 6, 3])
# 任意の部分集合を事象Xと呼ぶ
Out[6]:
In [26]:
"""
確率P(X) = 事象Xが起こる場合の数 / 起こりうる全ての場合の数
空事象Φ
P(7) = 0
余事象Ec
E = {1, 2, 3} のとき
Ec = {4, 5, 6}
積事象 A∩B キャップ
A = {1, 2, 3}
B = {1, 3, 4, 5} のとき
A∩B = {1, 3}
和事象 A∪B カップ
A∪B = {1,2,3,4,5}
"""
Out[26]:
In [37]:
# 統計的確率
"""サイコロをn回振る
根元事象が本当に1/6になっているか確認する
"""
# 1000回試行する
shikou = np.random.choice(dice_value, 1000)
# 1が出る回数を数える
bool_li = shikou == 1 # bool型のリストが得られる
kakuritu = bool_li.sum() / len(bool_li) # Trueの数 / 試行回数の数
print(kakuritu)
In [44]:
# 上記確率算出式から全根元事象について調べる
for i in range(1,7):
bool_li = shikou == i # bool型のリストが得られる
kakuritu = bool_li.sum() / len(bool_li) # Trueの数 / 試行回数の数
print("dice : ", i, ", probability:", kakuritu)
In [46]:
# 条件付き確率と乗法定理
"""
事象Aが生じた条件のもとで事象Bが生じる確率 = 「Aが与えられたもとでのBの条件付き確率」 と呼ぶ
P(B|A) = P(A∩B) / P(A)
変形させると乗法定理というものになる。
P(A∩B) = P(B|A) x P(A)
"""
Out[46]:
In [47]:
"""
例）サイコロを１回振って偶数が出た。その時 ４以上である確率は？
A = {2,4,6} # 偶数である条件
B = {4,5,6} # そのとき ４以上である場合の部分集合
条件付き確率を使うと
P(B|A) = P(A∩B) / P(A)
= P({4,6}) / P({2,4,6})
= (2/6) / (3/6)
= 2/3
"""
Out[47]:
In [48]:
"""
独立と従属
独立の場合はAという条件が付いても付かなくても同じ確率のはず
P(B|A) = P(B) ---(1)
また条件付き確率の式を使うと
P(B|A) = P(A∩B) / P(A) より
P(A∩B) = P(B|A) x P(A)
これに(1)を代入すると
P(A∩B) = P(B) x P(A) のはず
イコールでないとき 独立ではなく従属ということになる
例えば
左辺 = P(A∩B) = P({4,6}) = 2/6
右辺 = P(B) x P(A) =3/6 x 3/6
よって イコールにはならないのでAとBは従属関係であることがわかる。
"""
Out[48]:
In [49]:
"""
ベイズの定理
A:結果の事象
B:その原因の事象
P(B):Aが起きる前の確率（事前確率）
P(B|A): Aが起きたあとの事後Bの確率（事後確率）
P(A|B):Aが観測されたときにBが原因であるだろう確率（尤度）
P(B|A) = P(A|B) P(B) / [ P(A|B) P(B) + P(A|Bc) P(Bc) ]
"""
Out[49]:
In [50]:
"""練習問題"""
Out[50]:
In [51]:
"""確率変数と確率分布
用語：確率変数 実現値 離散確率変数 連続確率変数
分布関数 F(x)（累積確率分布関数）
密度関数 f(x)=F(x)/dx（確率密度関数）：分布関数の導関数
期待値 E(x)（平均）: Σ[x f(x)]
"""
Out[51]:
In [67]:
# 一様分布
# サイコロを１０００回振る
dice = np.array([1, 2, 3, 4, 5, 6])
choice = np.random.choice(dice, 1000)
# iが出る確率を計算する
kakuritu_li = []
for i in range(1,7):
bool_li = choice == i
kakuritu = bool_li.sum() / len(bool_li)
kakuritu_li.append(kakuritu)
# 確率を可視化する
plt.bar(dice, kakuritu_li)
plt.ylabel('probability')
Out[67]:
In [80]:
# ベルヌーイ分布
# コインを投げたとき tail = 0 , head = 1 とする
# 結果が[0,0,0,0,0,,1,1,1]のときの確率分布(headとtail)
coin_result = np.array([0, 0, 0, 0, 0, 1, 1, 1])
coin_unique = np.unique(coin_result) # np.array([0,1])と定義することと同じ
prob_li =[]
for i in coin_unique:
prob = sum(coin_result == i) / len(coin_result)
print(i, prob)
prob_li.append(prob)
print(prob_li)
# 可視化する
plt.bar(coin_unique, prob_li)
plt.xticks([0,1], ['tail', 'head'])
Out[80]:
In [94]:
# 二項分布 = 独立なベルヌーイ試行をn回繰り返したもの
# 二項分布をランダムに生成
x = np.random.binomial(10, 0.5, 10000)
plt.hist(x)
Out[94]:
In [95]:
# ポアソン分布 = 稀な事象が起きる確率
x = np.random.poisson(10, 10000)
plt.hist(x)
Out[95]:
In [98]:
# 正規分布
x = np.random.normal(5, 10, 10000) # 平均 標準偏差 サンプル数
plt.hist(x)
Out[98]:
In [102]:
# 対数正規分布
x = np.random.lognormal(5, 0.1, 10000)
plt.hist(x)
Out[102]:
In [120]:
# カーネル密度関数 = データから密度関数を推定する
import requests
import zipfile
from io import StringIO
import io
# データをurlからダウンロードする
url = 'http://archive.ics.uci.edu/ml/machine-learning-databases/00356/student.zip'
r = requests.get(url, stream=True)
z = zipfile.ZipFile(io.BytesIO(r.content))
z.extractall()
In [121]:
# データセットを読み込む
data_set = pd.read_csv('student-mat.csv', sep=';')
In [127]:
# カーネル密度関数
data_set.absences.plot(kind='kde', style='k--')
data_set.absences.hist(density=True)
plt.grid(True)
In [128]:
# 多次元確率分布
"""
用語：同時確率関数 周辺確率関数
条件付き確率関数 条件付き期待値
"""
Out[128]:
In [146]:
# ２次元の正規分布
import scipy.stats as st
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
# データの設定
x, y = np.mgrid[-100:100:2, -100:100:2]
pos = np.empty(x.shape+(2,))
pos[:, :, 0] = x
pos[:, :, 1] = y
In [159]:
# x y の平均が 0,0 共分散が 100,0 および 0,100
rv = multivariate_normal([0,0], [[100,0], [0,100]])
# 確率密度関数
z = rv.pdf(pos)
# 可視化
fig = plt.figure(dpi=100)
ax = Axes3D(fig)
ax.plot_wireframe(x,y,z)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('f(x,y)')
# 指数表記にする
ax.ticklabel_format(style='sci', axis='z', scilimits=(0,0))
In [1]:
# 推計統計学
"""
用語：大数の法則 中心極限定理
カイ２乗分布 スチューデントt分布 F分布
"""
Out[1]:
In [39]:
# 大数の法則
# サイコロを何回も振って平均値を出していくと期待値に近づいていく法則
def law():
dice = np.array([1, 2, 3, 4, 5, 6])
# 1000回試行する
result_li = [] # 出た目を格納するリスト
average_li = [] # 平均を格納するリスト
for i in range(10000):
result_li.append(np.random.choice(dice))
average = sum(result_li) / len(result_li) # 平均を算出する
average_li.append(average)
# 可視化
plt.plot(average_li)
plt.grid(True)
# ４回繰り返す
for shikou in range(4):
law()
In [41]:
# 別の方法 こっちのほうがリストのappendよりも高速
dice = np.array([1, 2, 3, 4, 5, 6])
shikou = np.arange(1,10001) # 1 から1000のリストを作成する。
for shikou_kaisuu in range(4):
choice = np.random.choice(dice, 10000) # dice から 1000回抽出する
# choiceを累積する
choice_cumsum = choice.cumsum()
# choice_cumsumを試行回数で除算して 平均を算出する
average = choice_cumsum / shikou
plt.plot(average)
plt.grid(True)
In [84]:
# 中心値 極限定理
"""
サイコロを振る回数が増えると標本平均が正規分布の形になる定理
"""
dice = np.array([1, 2, 3, 4, 5, 6])
average_li = []
def shikou(n=3):# 引数 サイコロを振る回数(１セット)
# サイコロをn回振る
choice = np.random.choice(dice, n)
# 平均を算出する
average = choice.sum() / n
#print(average)
average_li.append(average)
# １セットを何回も試行する
for cnt in range(10000):
shikou(3)
# 可視化 ヒストグラム
plt.figure(figsize=(10,5))
plt.subplot(1,3,1)
plt.hist(average_li,bins=20, range = (0, 8))
# shikou回数を変えてみる
for cnt in range(10000):
shikou(10)
plt.subplot(1,3,2)
plt.hist(average_li,bins=20, range = (0, 8))
# shikou回数を変えてみる
for cnt in range(10000):
shikou(100)
plt.subplot(1,3,3)
plt.hist(average_li,bins=20, range = (0, 8))
Out[84]:
In [105]:
# カイ２乗分布
"""
Zi : 標準正規分布（平均=0 σ=1）に従う確率変数
その２乗和を確率変数W W = Σ(Zi)^2 , i=1 to m とすると
これを自由度mのカイ事情分布に従う確率変数Wと呼ぶ
自由度nのカイ２乗分布
"""
def cisquare(m): # m = 自由度
x = np.random.chisquare(m, 1000)
return x
plt.hist(cisquare(5), color='blue')
plt.hist(cisquare(10), color='yellow')
plt.hist(cisquare(100), color='red')
Out[105]:
In [108]:
# スチューデントのt分布
"""
ZとWは独立な確率変数
Z:標準正規分布に従う
W:自由度mのカイ２乗分布に従う
とする
t分布 T=Z/sqrt(W/m)
"""
x = np.random.standard_t(5, 1000)
plt.hist(x)
Out[108]:
In [122]:
# スネディッカーのF分布
"""
W1:独立な確率変数, 自由度m1のカイ２乗分布に従う
W2:独立な確率変数, 自由度m2のカイ２乗分布に従う
F=(W1/m2)/(W2/m1)
とおいたとき、Fは自由度(m1,m2)のスネディッカーのF分布に従うという。
"""
plt.figure(figsize=(10,5))
#plt.subplot(1,2,1)
x = np.random.f(10, 10, 10000)
plt.hist(x, bins=10)
plt.grid(True)
#plt.subplot(1,2,2)
x = np.random.f(50, 50, 10000)
plt.hist(x, bins=10)
plt.grid(True)
In [123]:
#統計的推定
"""
用語：推定量 点推定
不偏性 一致性
"""
Out[123]:
In [124]:
"""
Xi:標本
標本平均: X_bar = (1/n)ΣXi
推定量：これを一般化すると確率変数の関数として記述できる
X_bar = T(X1,X2,...., Xn)
点推定：標本に基づいて母数を１点のパラメータとして言い当てること
不偏性 一致性：パラメータをより正確に推定できるかどうかの判断基準
推定量の期待値が母数θと一致するとき不偏性があるという・
E[T(X1,X2,Xn)] = θ
一致性：θの推定値E[T(X1,X2,Xn)]が観測個数nが大きくなるにつれてθに近づく性質
lim P[|T(X1,X2,Xn)-θ|>=ε] = 0
"""
Out[124]:
In [125]:
"""
区間推定
点推定では母数を１点で求めた
→区間推定はある程度の区間をもたせて母数を推定するもの
【説明】
平均がμ 分散が1の正規分布N(μ,1)から標本X1,X2,...,Xn が抽出されたとします。
この標本から母平均μを区間推定してみます。
標本平均X_bar は 平均μ分散1/nの正規分布N(μ,1/n)に従う ###なんで1/nなのか？
※標本のn数が多いほど分散(1/n)が小さくなるので、精度を上げるためにはn数はたくさん集めること。
α点：区間推定のしきい値？
Zα/2：正規分布のα点をこう書くらしい α/2はただの添字
P(-Zα/2 =< sqrt(n)(X_bar-μ) =< Zα/2) = 1-α
書き換えると
P(-Zα/2/sqrt(n) -X_bar =< (-μ) =< Zα/2/sqrt(n) -X_bar) = 1-α
P(+Zα/2/sqrt(n) +X_bar => (+μ) => -Zα/2/sqrt(n) +X_bar) = 1-α
よって
P(X_bar - Zα/2/sqrt(n) =< μ =< X_bar + Zα/2/sqrt(n)) = 1-α
信頼区間：上記範囲[X_bar - Zα/2/sqrt(n), X_bar + Zα/2/sqrt(n)] のこと
信頼係数：上記右辺 1-α のこと 信頼区間に母数（母平均）が入っている確率のこと
【以下はあまり使わないので必要なときに再度学ぶこととする。
最尤法
尤度関数
対数尤度関数
尤度方程式
尤度推定量
ベイズ法
"""
Out[125]:
In [2]:
"""
#統計的検定
#検定
#対立仮説
#棄却
#有意水準
#優位
#第一種の過誤
#第二種の過誤
#検出力
"""
Out[2]:
In [12]:
# データを読み込む
data_math = pd.read_csv('student-mat.csv', sep=';')
data_por = pd.read_csv('student-por.csv', sep=';')
# データをマージする
data_merge = pd.merge(data_math, data_por,
on=['school', 'sex', 'age', 'address', 'famsize', 'Pstatus', 'Medu', 'Fedu',
'Mjob', 'Fjob', 'reason', 'nursery', 'internet']
, suffixes=('_math', '_por'))
data_merge.head(3)
Out[12]:
In [17]:
# 平均を算出する
data_merge.keys()
print(data_merge.G1_math.mean())
print(data_merge.G1_por.mean())
In [21]:
#数学のほうが悪いように見えるが 本当にそうなのか？
#検定してみる
"""
帰無仮説H0：正しいか検討する仮説のこと。
（最終的に否定したい事柄をこれにする）
今回は母集団において平均に差がない、を帰無仮説としておく。
（否定して、差があると言いたい）
μmath = μpor
対立仮説H1：帰無仮説H0の否定
μmath ≠ μpor
ここで、帰無仮説H0 μmath = μpor が起こる確率が５％未満になるなら
この帰無仮説H0は棄却されるという。
そして、対立仮説H1が採択され、平均に差があるとみなされる。
このときの５％は有意水準と呼ばれる。α＝5% と書く
帰無仮説H0が正しい場合に、誤って帰無仮説H0を棄却してしまう確率P(reject|H0)のこと
p-value(p値)（有意確率）：
偶然 実際に反した数値が統計量として計算されてしまう確率。
例えばp値=0.002だった場合。帰無仮説H0が正しいなら今回起こったことは0.002の確率しか起きない現象ということ。
そんな低い確率が簡単に発生するとは考えにくいので、
そもそも帰無仮説H0は正しくない、という主張。
"""
# p値を算出する。
from scipy import stats
t, p = stats.ttest_rel(data_merge.G1_math, data_merge.G1_por)
print(p)
In [19]:
#第一種の過誤
"""
有意水準5%で帰無仮説H0が棄却されたとする。
しかし、もしかすると帰無仮説H0が正しかった場合もあり得る。
帰無仮説H0が正しかったにも関わらず棄却してしまうことを第一種の過誤という。
「正しいのに間違っていると言ってしまう」マチガイ
→偽陰性：α
あわてものマチガイ（疑い深く結論を急いでいるから）
信頼係数：（ 1-α ） 帰無仮説が正しいときに正しく採択できる確率を表す
"""
#第二種の過誤
"""
逆に、帰無仮説が
「間違っているのに正しいと言ってしまう」マチガイ
→第一種の過誤
偽陽性：β
ぼんやりもののマチガイ（見過ごしてしまっているから）
検出力：（１－β） 帰無仮説が間違っているときに正しく棄却できる確率を表す
"""
Out[19]:
In [41]:
# 練習問題 4-13
"""
data_math と data_por のキーG2について、平均に差があるか調べる。
帰無仮説H0:data_math.G2_math.mean() = data_math.G2_por.mearn()
有意確率p-valueを算出し、0.05以下であれば帰無仮説H0を棄却し平均に差があると判断する。
"""
from scipy import stats
# データを読み込む
data_math = pd.read_csv('student-mat.csv', sep=';')
data_por = pd.read_csv('student-por.csv', sep=';')
# データをマージする
data_merge = pd.merge(data_math, data_por,
on=['school', 'sex', 'age', 'address', 'famsize', 'Pstatus', 'Medu', 'Fedu',
'Mjob', 'Fjob', 'reason', 'nursery', 'internet']
, suffixes=('_math', '_por'))
data_merge.head(3)
#print(data_merge.keys())
# 有意確率p-valueを算出する
t, p = stats.ttest_rel(data_merge.G2_math, data_merge.G2_por)
print('data_merge.G2_math.mean() : ', data_merge.G2_math.mean())
print('data_merge.G2_por.mean() : ', data_merge.G2_por.mean())
print("p-value : ", p)
"""
「帰無仮説H0:平均値には差がない」有意確率p-value=4e-19であり、
有意差1%にも満たず極めて起こり得ないことである。
よって、「帰無仮説H0:平均値には差がない」を棄却し、
「平均値には差がある」と判断する。
"""
Out[41]:
In [42]:
# キーG3についても調べてみる
t, p = stats.ttest_rel(data_merge.G3_math, data_merge.G3_por)
print('data_merge.G3_math.mean() : ', data_merge.G3_math.mean())
print('data_merge.G3_por.mean() : ', data_merge.G3_por.mean())
print(p)
"""
「帰無仮説H0:平均値には差がない」有意確率p-value=5e-21であり、
有意差1%にも満たず極めて起こり得ないことである。
よって、「帰無仮説H0:平均値には差がない」を棄却し、
「平均値には差がある」と判断する。
"""
In [43]:
# キーabsencesについても調べてみる
t, p = stats.ttest_rel(data_merge.absences_math, data_merge.absences_por)
print('data_merge.absences_math.mean() : ',data_merge.absences_math.mean())
print('data_merge.absences_por.mean() : ',data_merge.absences_por.mean())
print(p)
"""
「帰無仮説H0:平均値には差がない」有意確率p-value=2e-6であり、
有意差1%にも満たず極めて起こり得ないことである。
よって、「帰無仮説H0:平均値には差がない」を棄却し、
「平均値には差がある」と判断する。
"""
In [56]:
# キーstudytimeについても調べてみる
t, p = stats.ttest_rel(data_merge.studytime_math, data_merge.studytime_por)
print('stats.ttest_rel(data_merge.studytime_math.mean() : ', data_merge.studytime_math.mean())
print('stats.ttest_rel(data_merge.studytime_por.mean() : ', data_merge.studytime_por.mean())
print(p)
"""
「帰無仮説H0:平均値には差がない」有意確率p-value=0.56であり、
有意差56%にも満たず極めて起こり得ないことである。
よって、「帰無仮説H0:平均値には差がない」を棄却することはできないので、
「平均値には差がない」と判断する。
"""