データサイエンティスト養成講座4章の演習

引き続きデータサイエンティスト養成講座を読む。

東京大学のデータサイエンティスト育成講座

作者:塚本邦尊,山田典一,大澤文孝
マイナビ出版

Amazon

4章では、数理統計学の初歩的な部分とその実装について学ぶ。数理統計について最低限の知識がないと理解がしづらいと思う。自分は数理統計の初歩的な部分を勉強するためにクロスセクションシリーズの「数理統計学の基礎」を読んだ。

数理統計学の基礎 (クロスセクショナル統計シリーズ 1)

作者:尾畑伸明
共立出版

Amazon

この本は数理統計学の入門書の中では比較的、厳密に書かれていると思う。

演習問題4-1

コイン投げの試行を1000回実行し、成功と失敗の確率を示せ

coin_data = np.array([0, 1])
steps = 1000
coin_rolls = np.random.choice(coin_data, steps)
for i in range(0, 2):
  p = len(coin_rolls[coin_rolls==i])/steps
  if(i == 0):
    print(f'表が出る確率は{p}')
  else:
    print(f'裏が出る確率は{p}')

演習問題4-2

1000本のくじの中に100本の当たりがある。AとBが順番にくじを引くときにともにあたりが出る確率を求めよ(くじは元に戻さない)

A、Bが当たりを引く確率をそれぞれ $P(A)$ 、 $P(B)$ と書く。 $P(A \cap B ) = P(A) \, P(B|A) = \frac{100}{1000} \times \frac{99}{999} = \frac{11}{1110}$

p.108のカーネル密度推定は他の本で勉強をする。

練習問題4-4

N(0, 1)からn=100の標本抽出を10000回繰り返して標本平均のヒストグラムを図示

sample_mean = np.array([])
for _ in range(0, 10000):
  #標本抽出して、標本平均を求める
  x = np.random.normal(0, 1, 100)
  sample_mean = np.append(sample_mean, x)

plt.hist(sample_mean)
plt.grid(True)

練習問題4-5 練習問題4-4の対数正規分布の場合

sample_mean = np.array([])
for _ in range(0, 10000):
  #標本抽出して、標本平均を求める
  x = np.random.lognormal(0, 1, 100)
  sample_mean = np.append(sample_mean, x)
#標本平均はほとんど0から15の範囲に分布しているので、ヒストグラムの範囲を調整する
plt.hist(sample_mean,range=(0, 15))
plt.grid(True)

練習問題4-6 student_data_mathの一期目の成績G1のヒストグラムとカーネル密度推定

import pandas as pd
student_data_math = pd.read_csv('student-mat.csv', sep=';')

#カーネル密度推定
student_data_math.G1.plot(kind='kde', style='k--')

#ヒストグラム
student_data_math.G1.hist(density=True)
plt.grid(True)

カーネル密度推定とヒストグラム

p.111の二次元正規分布の図示

ライブラリの機能とかをまだ理解していないので、実装の意味が分からないところが結構ある pos = np.empty(x.shape + (2, ))、pos[:, :, 0] = x pos[:, :, 1] = yとかよくわからん

#2次元正規分布
import scipy.stats as st
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D

x, y = np.mgrid[10:100:2, 10:100:2]

#print(x)
#print(y)

pos = np.empty(x.shape + (2, ))
print(pos.shape)
pos[:, :, 0] = x
pos[:, :, 1] = y
print(pos[0])

rv = multivariate_normal([50,50], [[100, 0], [0, 100]])
z = rv.pdf(pos)
print(z)

fig = plt.figure(dpi=100)

ax = Axes3D(fig)
ax.plot_wireframe(x, y, z)
#x,y,zラベルの設定
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

ax.ticklabel_format(style='sci', axis='z', scilimits=(0, 0))

大数の法則

大数の法則はiidな確率変数 $X _ 1, X _ 2, \cdots X _ n$ について $\displaystyle \sum _ {k=1}^{n} \frac{x _ {k}}{n}$ が $n \rightarrow \infty$ としたときにpに確率収束することを示している。

#大数の法則
calc_times = 1000
#サイコロ
sample_array = np.array([1,2,3,4,5,6])
number_cnt =  np.arange(1, calc_times+1)

#4つのパスを作成
for i in range(4):
  #1から6から一つ抽出する操作をcalc_times回繰り返してその累積和をとる
  p = np.random.choice(sample_array, calc_times).cumsum()
  #print(p)
  plt.plot(p/number_cnt)
  plt.grid(True)

中心極限定理

#中心極限定理
def function_central_theory(N):

  sample_array = np.array([1, 2, 3, 4, 5, 6])
  number_cnt = np.arange(1, N+1) * 1.0
  
  mean_array = np.array([])

  for _ in range(1000):
    cum_variables = np.random.choice(sample_array, N).cumsum()*1.0
    print(cum_variables[N-1])
    mean_array = np.append(mean_array, cum_variables[N-1]/N)
  
  plt.hist(mean_array)
  plt.grid(True)

function_central_theory(100)

練習問題4-10は省略(どの教科書にも書いてあって、明らかなので)

練習問題4-11

Pythonで最尤推定の実装をする際に参考にする

練習問題4-14

本と同様にp値を計算するコードを書いたけど、計算できてない

#p値を計算
from scipy import stats
#G2での数学とポルトガル語の成績の平均の差の検定
t_2, p_2 = stats.ttest_rel(student_data_merge.G2_math, student_data_merge.G2_por)
print("p値は, ", p_2)

ながれとよどみ

自分が解いた問題の解法など