Warning: Undefined variable $position in /home/pystyles/pystyle.info/public_html/wp/wp-content/themes/lionblog/functions.php on line 4897
  1. HOME
  2. Python
  3. scipy - scipy.stats の使い方まとめ

scipy – scipy.stats の使い方まとめ

scipy – scipy.stats の使い方まとめ
目次

概要

確率分布を扱う scipy のモジュールである scipy.stats の使い方について解説します。

Advertisement

確率変数

離散型及び連続型確率変数を表すクラスが提供されています。

scipy.stats.rv_discrete

scipy.stats.rv_discrete は離散型確率変数を表すクラスです。

In [1]: 
import numpy as np
from matplotlib import pyplot as plt
import scipy.stats as stats

# ポアソン分布に従う確率変数を作成する。
X = stats.poisson(2)

サンプリング

scipy.stats.rv_discrete.rvs() で確率変数から観測値を取得します。

In [2]: 
print(X.rvs(size=20))
[4 1 0 2 1 1 3 1 2 0 2 1 1 3 3 2 3 1 0 1]

確率質量関数

scipy.stats.rv_discrete.pmf() は、確率質量関数の値 $f_X(x)$ を返します。

In [3]: 
fig, ax = plt.subplots()

x = np.arange(20)
y = X.pmf(x)

ax.stem(x, y, use_line_collection=True)
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

scipy.stats.rv_discrete.logpmf() は、確率質量関数の対数をとった値 $\log f_X(x)$ を返します。

In [4]: 
fig, ax = plt.subplots()

x = np.arange(20)
y = X.logpmf(x)

ax.stem(x, y, use_line_collection=True)
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

累積分布関数

scipy.stats.rv_discrete.cdf() は、累積分布関数 $F_X(x) = P(X \le x)$ の値を返します。

In [5]: 
fig, ax = plt.subplots()

x = np.linspace(0, 20, 100)
y = X.cdf(x)

ax.plot(x, y)
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

scipy.stats.rv_discrete.logcdf() は、累積分布関数の対数をとった値 $\log F_X(x)$ を返します。

In [6]: 
fig, ax = plt.subplots()

x = np.linspace(0, 20, 100)
y = X.logcdf(x)

ax.plot(x, y)
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

生存関数

scipy.stats.rv_discrete.sf() は、生存関数 $S_X(x) = P(X > x) = 1 – F_X(x)$ の値を返します。

In [7]: 
fig, ax = plt.subplots()

x1 = np.linspace(0, 20, 100)
y1 = X.cdf(x1)

x2 = np.linspace(0, 20, 100)
y2 = X.sf(x2)

ax.plot(x1, y1, label="cdf")
ax.plot(x2, y2, label="sf")
ax.legend()
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

scipy.stats.rv_discrete.logsf() は、生存関数の対数をとった値 $\log S_X(x)$ を返します。

In [8]: 
fig, ax = plt.subplots()

x1 = np.linspace(0, 20, 100)
y1 = X.logcdf(x1)

x2 = np.linspace(0, 20, 100)
y2 = X.logsf(x2)

ax.plot(x1, y1, label="logcdf")
ax.plot(x2, y2, label="logsf")
ax.legend()
ax.grid()
ax.set_xticks(np.arange(0, 21, 2))

plt.show()

分位関数 (累積分布関数の逆関数)

scipy.stats.rv_discrete.ppf() は分位関数 $Q(p) = \inf \{ x \in \R \mid p \le F(x) \}$ の値を返します。

In [9]: 
fig, ax = plt.subplots()

p = np.linspace(0, 1, 100)
x = X.ppf(p)

ax.plot(p, x)
ax.grid()

plt.show()

生存関数の逆関数

scipy.stats.rv_discrete.isf() は生存関数の逆関数の値 $Q(p) = \inf \{ x \in \R \mid p \le S_X(x) \}$ を返します。

In [10]: 
fig, ax = plt.subplots()

p = np.linspace(0, 1, 100)
x = X.isf(p)

ax.plot(p, x)
ax.grid()

plt.show()

モーメント

scipy.stats.rv_discrete.moment() は $n$ 次のモーメントを計算します。

In [11]: 
for i in range(1, 5):
    m = X.moment(i)
    print(fr"E[X^{i}] = {m:.2f}")
E[X^1] = 2.00
E[X^2] = 6.00
E[X^3] = 22.00
E[X^4] = 94.00

平均情報量 / エントロピー

scipy.stats.rv_discrete.entropy() で平均情報量を計算します。

In [12]: 
print(X.entropy())
1.7048826439329838

統計量

scipy.stats.rv_discrete.stats() で次の統計量を計算できます。計算対象の統計量は文字列で指定します。

  • “m”: 期待値
  • “v”: 分散
  • “s”: 歪度 (skewness)
  • “k”: 尖度 (kutosis)
In [13]: 
mean, var, skew, kutosis = X.stats("mvsk")
print(f"mean = {mean:.2f}")
print(f"var = {var:.2f}")
print(f"skew = {skew:.2f}")
print(f"kutosis = {kutosis:.2f}")
mean = 2.00
var = 2.00
skew = 0.71
kutosis = 0.50
In [14]: 
print(f"median = {X.median():.2f}")
print(f"mean = {X.mean():.2f}")
print(f"var = {X.var():.2f}")
print(f"std = {X.std():.2f}")
median = 2.00
mean = 2.00
var = 2.00
std = 1.41

確率変数を引数にとる関数 $f$ の期待値 $E[f(X)]$ を計算する場合は scipy.stats.rv_discrete.expect() を使用します。

In [15]: 
def f(x):
    return x ** 2 + 1


print(X.expect(f))
6.999999999999999

信頼区間

scipy.stats.rv_discrete.interval() は、$\alpha$ が与えられたとき、$P(median[X] – t \le X \le median[x] + t) = \alpha$ を満たす区間 $[median[X] – t, median[x] + t]$ を計算します。

In [16]: 
a, b = X.interval(0.95)
print(f"P[{a}, {b}]")
P[0.0, 5.0]

台

scipy.stats.rv_discrete.support() は、確率質量関数の台 $\text{supp}(pmf) = \{x \in \R \mid f_X(x) > x \}$ を返します。

In [17]: 
print(X.support())
(0, inf)

scipy.stats.rv_continuous

scipy.stats.rv_continuous は連続型確率変数を表すクラスです。

In [18]: 
# 平均0、分散1に従う確率変数を作成する。
X = stats.norm(0, 1)

loc, scale パラメータ

すべての連続型確率変数の分布は、位置やスケールを調整するためのパラメータ loc, scale をとります。 調整前の確率変数が $X$ としたとき、調整後の確率変数は以下のように計算されます。

$X’ = \frac{X – loc}{scale}$

In [19]: 
fig, ax = plt.subplots()

x = np.linspace(-10, 10, 100)
y1 = stats.norm.pdf(x)
y2 = stats.norm.pdf(x, loc=2, scale=2)

ax.plot(x, y1)
ax.plot(x, y2)
ax.grid()

plt.show()

サンプリング

scipy.stats.rv_continuous.rvs() で確率変数から観測値を取得します。

In [20]: 
print(X.rvs(size=20))
[ 0.52587614 -1.07375603 -1.24883296 -0.59328601  0.17404647 -0.30302527
 -1.12587053  0.47764943 -0.81734939  0.16761539 -0.09701077  1.01066471
  0.58933414  1.44496678 -0.34886864 -2.59756075  0.07253381  1.07815998
 -0.56592292  0.26258007]

確率密度関数

scipy.stats.rv_continuous.pdf() は、確率密度関数の値 $f_X(x)$ を返します。

In [21]: 
fig, ax = plt.subplots()

x = np.linspace(-5, 5, 100)
y = X.pdf(x)

ax.plot(x, y)
ax.grid()

plt.show()

scipy.stats.rv_continuous.logpdf() は、確率密度関数の対数をとった値 $\log f_X(x)$ を返します。

In [22]: 
fig, ax = plt.subplots()

x = np.linspace(-5, 5, 100)
y = X.logpdf(x)

ax.plot(x, y)
ax.grid()

plt.show()

累積分布関数

scipy.stats.rv_continuous.cdf() は、累積分布関数 $F_X(x) = P(X \le x)$ の値を返します。

In [23]: 
fig, ax = plt.subplots()

x = np.linspace(-5, 5, 100)
y = X.cdf(x)

ax.plot(x, y)
ax.grid()

plt.show()

scipy.stats.rv_continuous.logcdf() は、累積分布関数の対数をとった値 $\log F_X(x)$ を返します。

In [24]: 
fig, ax = plt.subplots()

x = np.linspace(-5, 5, 100)
y = X.logcdf(x)

ax.plot(x, y)
ax.grid()

plt.show()

生存関数

scipy.stats.rv_continuous.sf() は、生存関数 $S_X(x) = P(X > x) = 1 – F_X(x)$ の値を返します。

In [25]: 
fig, ax = plt.subplots()

x1 = np.linspace(-5, 5, 100)
y1 = X.cdf(x1)

x2 = np.linspace(-5, 5, 100)
y2 = X.sf(x2)

ax.plot(x1, y1, label="cdf")
ax.plot(x2, y2, label="sf")
ax.legend()
ax.grid()

plt.show()

scipy.stats.rv_continuous.logsf() は、生存関数の対数をとった値 $\log S_X(x)$ を返します。

In [26]: 
fig, ax = plt.subplots()

x1 = np.linspace(-5, 5, 100)
y1 = X.logcdf(x1)

x2 = np.linspace(-5, 5, 100)
y2 = X.logsf(x2)

ax.plot(x1, y1, label="logcdf")
ax.plot(x2, y2, label="logsf")
ax.legend()
ax.grid()

plt.show()

分位関数 (累積分布関数の逆関数)

scipy.stats.rv_continuous.ppf() は分位関数 $Q(p) = \inf \{ x \in \R \mid p \le F(x) \}$ の値を返します。

In [27]: 
fig, ax = plt.subplots()

p = np.linspace(0, 1, 100)
x = X.ppf(p)

ax.plot(p, x)
ax.grid()

plt.show()

生存関数の逆関数

scipy.stats.rv_continuous.isf() は生存関数の逆関数の値 $Q(p) = \inf \{ x \in \R \mid p \le S_X(x) \}$ を返します。

In [28]: 
fig, ax = plt.subplots()

p = np.linspace(0, 1, 100)
x = X.isf(p)

ax.plot(p, x)
ax.grid()

plt.show()

モーメント

scipy.stats.rv_continuous.moment() は $n$ 次のモーメント $E[X^n]$ を計算します。

In [29]: 
for i in range(1, 5):
    m = X.moment(i)
    print(fr"E[X^{i}] = {m:.2f}")
E[X^1] = 0.00
E[X^2] = 1.00
E[X^3] = 0.00
E[X^4] = 3.00

平均情報量 / エントロピー

scipy.stats.rv_continuous.entropy() で平均情報量を計算します。

In [30]: 
print(X.entropy())
1.4189385332046727

統計量

scipy.stats.rv_continuous.stats() で次の統計量を計算できます。計算対象の統計量は文字列で指定します。

  • “m”: 期待値
  • “v”: 分散
  • “s”: 歪度 (skewness)
  • “k”: 尖度 (kutosis)
In [31]: 
mean, var, skew, kutosis = X.stats("mvsk")
print(f"mean = {mean:.2f}")
print(f"var = {var:.2f}")
print(f"skew = {skew:.2f}")
print(f"kutosis = {kutosis:.2f}")
mean = 0.00
var = 1.00
skew = 0.00
kutosis = 0.00
In [32]: 
print(f"median = {X.median():.2f}")
print(f"mean = {X.mean():.2f}")
print(f"var = {X.var():.2f}")
print(f"std = {X.std():.2f}")
median = 0.00
mean = 0.00
var = 1.00
std = 1.00

確率変数を引数にとる関数 $f$ の期待値 $E[f(X)]$ を計算する場合は scipy.stats.rv_continuous.expect() を使用します。

In [33]: 
def f(x):
    return x ** 2 + 1


print(X.expect(f))
2.000000000000001

信頼区間

scipy.stats.rv_continuous.interval() は、$\alpha$ が与えられたとき、$P(median[X] – t \le X \le median[x] + t) = \alpha$ を満たす区間 $[median[X] – t, median[x] + t]$ を返します。

In [34]: 
a, b = X.interval(0.95)
print(f"P[{a}, {b}]")

fig, ax = plt.subplots()

x1 = np.linspace(-5, 5, 100)
y1 = X.pdf(x1)

x2 = np.linspace(a, b, 100)
y2 = X.pdf(x2)

ax.plot(x1, y1)
ax.fill_between(x2, np.zeros_like(y2), y2, alpha=0.3)
ax.grid()

plt.show()
P[-1.959963984540054, 1.959963984540054]

台

scipy.stats.rv_continuous.support() は、確率質量関数の台 $\text{supp}(pdf) = \{x \in \R \mid f_X(x) > 0 \}$ を返します。

In [35]: 
print(X.support())
(-inf, inf)

最尤推定

scipy.stats.rv_continuous.fit() でその分布の最尤推定を行います。

In [36]: 
data = X.rvs(size=100)

mean, std = stats.norm.fit(data)
print(f"mean: {mean}, std: {std}")
mean: 0.18434298655426873, std: 1.038091568649953

scipy.stats.rv_continuous.fit_loc_scale() は loc, scale を推定します。

In [37]: 
data = stats.uniform.rvs(loc=5, scale=1, size=1000)

loc, scale = stats.uniform.fit_loc_scale(data)
print(f"loc: {loc}, scale: {scale}")
loc: 4.986848789248153, scale: 1.0205864431590375

対数尤度関数

scipy.stats.rv_continuous.nnlf() で対数尤度関数の負をとった値を返します。 theta には、分布を表すパラメータ (loc、scale を含む) を theta = (*args, loc, scale) とタプルで渡します。

In [38]: 
loc, scale = 0, 1
data = stats.norm.rvs(loc=loc, scale=scale, size=100)

nnl = stats.norm.nnlf((loc, scale), data)
print(nnl)
144.17743898141782