Pythonを用いた統計学の手法のプログラミング 〜NumPy、Matplotlib etc〜

数式だけの解説ではわかりにくい場合もあると思われるので、統計学の手法や関連する概念をPythonのプログラミングで表現します。

基本的な公式

平均

$$
\large
\begin{align}
\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
print(np.mean(x))

・実行結果

> print(np.mean(x))
5.0

分散

$V[X]=E[(X-E[X])^2]$を元に計算

$$
\large
\begin{align}
S^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i-\bar{x})^2
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
print(np.sum((x-np.mean(x))**2)/(x.shape[0]))

・実行結果

> print(np.sum((x-np.mean(x))**2)/(x.shape[0]))
2.0

$V[X]=E[X^2]-E[X]^2$を元に計算

$$
\large
\begin{align}
S^2 = \frac{1}{n} \sum_{i=1}^{n} x_i^2 – \bar{x}^2
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
print(np.sum(x**2)/(x.shape[0])-np.mean(x)**2)

・実行結果

> print(np.sum(x**2)/(x.shape[0])-np.mean(x)**2)
2.0

ここで用いた数式の導出は下記で行なった。
https://www.hello-statisticians.com/explain-terms-cat/expectation-variance-covariance.html

共分散

$Cov[X,Y]=E[(X-E[X])(Y-E[Y])]$を元に計算

$$
\large
\begin{align}
S_{XY} = \frac{1}{n} \sum_{i=1}^{n} (x_i-\bar{x})(y_i-\bar{y})
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
y = np.array([3.1, 7.1, 3.9, 5.9, 5.])

print(np.dot(x-np.mean(x),y-np.mean(y))/(x.shape[0]))
print(np.sum((x-np.mean(x))*(y-np.mean(y)))/(x.shape[0]))

・実行結果

> print(np.dot(x-np.mean(x),y-np.mean(y))/(x.shape[0]))
2.0
> print(np.sum((x-np.mean(x))*(y-np.mean(y)))/(x.shape[0]))
2.0

np.dotを用いて内積を計算する方法と、ベクトルの要素ごとの積にnp.sumを適用する方法の主に二通りの実装ができる。

$Cov[X,Y]=E[XY]-E[X]E[Y]$を元に計算

$$
\large
\begin{align}
S^2 = \frac{1}{n} \sum_{i=1}^{n} x_iy_i – \bar{x}\bar{y}
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
y = np.array([3.1, 7.1, 3.9, 5.9, 5.])

print(np.sum(x*y)/(x.shape[0])-np.mean(x)*np.mean(y))
print(np.dot(x,y)/(x.shape[0])-np.mean(x)*np.mean(y))

・実行結果

> print(np.sum(x*y)/(x.shape[0])-np.mean(x)*np.mean(y))
2.0
> print(np.dot(x,y)/(x.shape[0])-np.mean(x)*np.mean(y))
2.0

ここで用いた数式の導出は下記で行なった。
https://www.hello-statisticians.com/explain-terms-cat/expectation-variance-covariance.html

共分散行列

下記で作成を行った。
https://www.hello-statisticians.com/explain-books-cat/stat_workbook/stat_workbook_ch22.html#221

相関係数

$$
\large
\begin{align}
r_{XY} = \frac{\sum_{i=1}^{n} (x_i-\bar{x})(y_i-\bar{y})}{\sqrt{\sum_{i=1}^{n} (x_i-\bar{x})^2} \sqrt{\sum_{i=1}^{n} (y_i-\bar{y})^2}}
\end{align}
$$

上記の数式は、下記のようにPythonで表すことができる。

import numpy as np

x = np.array([3., 7., 4., 6., 5.])
y = np.array([3.1, 7.1, 3.9, 5.9, 5.])

r = np.dot(x-np.mean(x),y-np.mean(y))/np.sqrt(np.dot(x-np.mean(x),x-np.mean(x))*np.dot(y-np.mean(y),y-np.mean(y)))

print(r)

・実行結果

> print(r)
0.99800598007

相関行列

下記で作成を行った。
https://www.hello-statisticians.com/explain-books-cat/stat_workbook/stat_workbook_ch22.html#221

確率分布の描画

離散分布

ポアソン分布

$$
\large
\begin{align}
f(x|\lambda) = \frac{\lambda^{x} e^{-\lambda}}{x!}, \quad x \geq 0
\end{align}
$$

上記で表したポアソン分布$Po(\lambda)$の確率密度関数のグラフは下記を実行することで描画を行うことができる。

import math
import matplotlib.pyplot as plt
import numpy as np

def factorial(vec):
    res = np.zeros(vec.shape[0])
    for i in range(vec.shape[0]):
        res[i] = math.factorial(vec[i])
    return res

lamb1, lamb2, lamb3 = 3., 5., 7.

x = np.arange(0, 16, 1)
y1 = lamb1**x * np.e**(-lamb1) / factorial(x)
y2 = lamb2**x * np.e**(-lamb2) / factorial(x)
y3 = lamb3**x * np.e**(-lamb3) / factorial(x)

plt.plot(x,y1,color="blue")
plt.plot(x,y2,color="green")
plt.plot(x,y3,color="red")
plt.show()

・実行結果

二項分布の極限① ポアソン小数の法則とポアソン分布

「二項分布の極限のプログラミング」の「ポアソンの小数の法則」で取り扱いを行った。

二項分布の極限② 中心極限定理と正規分布

「二項分布の極限のプログラミング」の「中心極限定理」で取り扱いを行った。

連続分布

正規分布

$$
\large
\begin{align}
f(x|\mu,\sigma) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left[ -\frac{(x-\mu)^2}{2 \sigma^2} \right]
\end{align}
$$

上記で表した正規分布$N(\mu,\sigma)$の確率密度関数のグラフは下記を実行することで描画を行うことができる。

import numpy as np
import matplotlib.pyplot as plt

mu, sigma = 0, 1

x = np.arange(-5., 5.01, 0.01)
y = np.exp(-(x-mu)**2/2*sigma**2)/np.sqrt(2*np.pi*sigma**2)

plt.plot(x,y)
plt.show()

・実行結果

指数分布

$$
\large
\begin{align}
f(x|\lambda) = \lambda e^{-\lambda x}, \quad x \geq 0
\end{align}
$$

上記で表した正規分布$Exp(\lambda)$の確率密度関数のグラフは下記を実行することで描画を行うことができる。

import numpy as np
import matplotlib.pyplot as plt

lamb = 1

x = np.arange(0., 5.01, 0.01)
y = lamb * np.e**(-lamb*x)

plt.plot(x,y)
plt.show()

・実行結果

乱数生成

乗算合同法

「自然科学の統計学」の章末課題11.2の解答で取り扱った。

M系列

「自然科学の統計学」の章末課題11.3の解答で取り扱った。

メルセンヌツイスタ法

「自然科学の統計学」の章末課題11.1の解答で取り扱った。

中心極限定理を用いた正規乱数の生成

「自然科学の統計学」の章末課題11.5の解答で取り扱った。

ボックス・ミュラー法

確率過程・時系列解析

ブラウン運動

「確率過程のプログラミング」の「ブラウン運動」で取り扱いを行った。

ポアソン過程

「確率過程のプログラミング」の「ポアソン過程」で取り扱いを行った。

ARMA過程・共分散定常

下記で取り扱った。
https://www.hello-statisticians.com/python/stat_program3.html

ホワイトノイズ

コレログラム

統計的推定

統計数値表

「Pythonを用いた統計的推測のプログラミング」の「統計数値表とSciPyを用いた代用」で取り扱った。

区間推定

・分散既知の場合の母平均の区間推定

仮説検定

・分散既知の場合の母平均の検定

ノンパラメトリック法

順位和検定

「Pythonを用いたノンパラメトリック法のプログラミング」の「順位和検定」で取り扱った。

並べ替え検定

「統計学実践ワークブック」の「Ch.13 ノンパラメトリック法」の解答例で取り扱った。

符号付き順位検定

符号検定

「Pythonを用いたノンパラメトリック法のプログラミング」の「符号検定」で取り扱った。

クラスカル・ウォリス検定

「統計学実践ワークブック」の「Ch.13 ノンパラメトリック法」の解答例で取り扱った。

回帰

線形回帰

単回帰のパラメータ推定

「Pythonを用いた回帰のプログラミング」の「単回帰のパラメータ推定」で取り扱った。

重回帰のパラメータ推定

「Pythonを用いた回帰のプログラミング」の「重回帰のパラメータ推定」で取り扱った。

質的回帰

ロジスティック回帰のパラメータ推定

「Pythonを用いた回帰のプログラミング」の「ロジスティック回帰」で取り扱った。

重回帰のパラメータ推定

「Pythonを用いた回帰のプログラミング」の「プロビット回帰」で取り扱った。

離散値の回帰

ポアソン回帰

「Pythonを用いた回帰のプログラミング」の「ポアソン回帰」で取り扱った。

ベイズ法

共役事前分布

二項分布-ベータ分布

「Pythonを用いたベイズ法のプログラミング」の「二項分布とベータ分布」で取り扱った。

ポアソン分布-ガンマ分布

正規分布-正規分布

MCMC

メトロポリス法

「Pythonを用いたMCMCのプログラミング」の「メトロポリス法」で取り扱った。

自動微分

自動微分の基本処理

「Pythonを用いた自動微分のプログラミング」の「自動微分の基本処理」で取り扱った。

自動微分のモジュール化

「Pythonを用いた自動微分のプログラミング」の「自動微分のモジュール化」で取り扱った。

注意しておくとよい事項

変数の型がint型の際の割り算の結果

int型をint型で割る際は、結果が小数ではなく余りの切り捨ての形式になることは注意しておくと良い。

import numpy as np

x = np.array([1, 2, 3, 1, 2])

print("np.mean(x): {}".format(np.mean(x)))
print("np.sum(x)/5: {}".format(np.sum(x)/5))
print("np.sum(x)/5.: {}".format(np.sum(x)/5.))

・実行結果

> print("np.mean(x): {}".format(np.mean(x)))
np.mean(x): 1.8
> print("np.sum(x)/5: {}".format(np.sum(x)/5))
np.sum(x)/5: 1
> print("np.sum(x)/5.: {}".format(np.sum(x)/5.))
np.sum(x)/5.: 1.8