pythonでF分布の計算をする


F分布について、python のライブラリでの計算方法を調べてみました。
結果を記載します。


F分布について

そもそもポアソン分布とは何かの理解が曖昧なので、書籍、Webで調べた結果をざっくり記載します。
* 確率変数 U と V が互いに独立で、U は 自由度 N1 の カイ二乗分布に従い、 V は 自由度N2 の カイ二乗分布に従う場合、フィッシャーの分散比の取り得る分布を F分布という。
* F分布はF検定で使用される。 F検定では、データXおよびデータYで与えられた2群間における分散が等しいかどうかを判断できる。
* F検定は、等分散性の検定、重回帰分析のパラメータの検定 に使用される。

F検定を実施するために必要なものが、F分布だと個人的には理解しました。


F分布の計算ができるライブラリについて

  • np.random.f
  • scipy.stats.f

各ライブラリを使ってみる

np.random.f

numpy.random.f を使うと、F分布に従う、擬似乱数を生成することができます。
numpy.random.f(dfnum, dfden, size=None) の引数 dfnum は分子の自由度、 dfden は分母の自由度、size は戻り値の要素数を指定可能です。

# 自動度(3, 10) のF分布 乱数を生成    
import numpy as np
s = np.random.f(3., 10., 100)
# ヒストグラム描画 
%matplotlib inline
import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(s, 14, normed=True)
print("上側確率 5%:", sorted(s)[-5])
print("上側確率 1%:", sorted(s)[-1])
plt.show()


# 自動度(3, 10) のF分布 乱数を生成    
import numpy as np
s = np.random.f(3., 10., 1000)
# ヒストグラム描画 
%matplotlib inline
import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(s, 14, normed=True)
print("上側確率 5%:", sorted(s)[-50])
print("上側確率 1%:", sorted(s)[-10])
plt.show()

# 自動度(3, 10) のF分布 乱数を生成    
import numpy as np
s = np.random.f(3., 10., 10000)
# ヒストグラム描画 
%matplotlib inline
import matplotlib.pyplot as plt
count, bins, ignored = plt.hist(s, 14, normed=True)
print("上側確率 5%:", sorted(s)[-500])
print("上側確率 1%:", sorted(s)[-100])
plt.show()
上側確率 5%: 3.38566314789
上側確率 1%: 10.7104924276

png

上側確率 5%: 4.00906109216
上側確率 1%: 7.54627385611

png

上側確率 5%: 3.79496866783
上側確率 1%: 6.3994598002

png

F分布表 の、分母 10 、分子 3 の 5%点は、3.708 で、
F分布表 の、分母 10 、分子 3 の 1%点は、6.552 です。
結果を見る限り、うまくF分布の乱数の生成ができているように思います。

scipy.stats.f

scipy.stats.f を使い、F分布のグラフの描画、p値を求めることができます。

  • グラフを描画する
from scipy.stats import f
import matplotlib.pyplot as plt

%matplotlib inline
fig, ax = plt.subplots(1, 1)

# 分子、分母の自由度
dfn, dfd = 3, 10

# mements は 文字列 mvsk を指定すると、 m > 平均 v > 分散 >  s > 歪度(わいど)  k > 尖度(せんど) が返る    
# このあたりは、scipy.stats.poisson と使い方は同じ
mean, var, skew, kurt = f.stats(dfn, dfd, moments='mvsk')
print("平均", mean)
print("分散", var)
print("歪度", skew)
print("尖度", kurt)

# F分布を描画
x = np.linspace(f.ppf(0.01, dfn, dfd),f.ppf(0.99, dfn, dfd), 100)
ax.plot(x, f.pdf(x, dfn, dfd), 'r-', lw=5, alpha=0.6, label='f pdf')

# パラメータを固定?した F分布を描画
rv = f(dfn, dfd)
ax.plot(x, rv.pdf(x), 'k-', lw=2, label='frozen pdf')

# 2つの配列が近似的に同じか比較する
vals = f.ppf([0.001, 0.5, 0.999], dfn, dfd)
result = np.allclose([0.001, 0.5, 0.999], f.cdf(vals, dfn, dfd))
print("配列は近似的に同様か:", result)
平均 1.25
分散 1.9097222222222223
歪度 4.2211588240886915
尖度 59.45454545454547
配列は近似的に同様か: True

png

  • p値の計算をする
# scipy.statsから、 poissonをインポートします
from scipy.stats import f

# 分子、分母の自由度
dfn, dfd = 3, 10

# 5%点の値(上側なので、95% を指定する)
five = f.ppf(0.95, dfn, dfd)
print('上側確率 5%:', five)

# 1%点の値(上側なので、99% を指定する)
one = f.ppf(0.99, dfn, dfd)
print('上側確率 1%:', one)
上側確率 1%: 6.55231255752
上側確率 5%: 3.70826481905

参考

以下、参考になりました。

以上です。

コメント