【PythonのkstestでKS検定】出荷数が少なくても正規分布になることを検証する
観測データがある確率分布に従っているかどうかを検定するのにKS検定は用いられますが、中でも正規分布への適合性の判定によく用いられます。
安全在庫理論でも出荷データが正規分布に従うかどうかは大きな関心事です。
特に日々の出荷数が少ないC商品は、そのままでは正規分布にならない可能性が高いでしょう。
しかし、だからといって安全在庫理論が使えないということはありません。
日々の出荷データは正規分布にならなくても、週でまとめれば正規分布になることがほとんどです。
仮に駄目でも、月でまとめれば正規分布になるでしょう。
これをKS検定で調べてみましょう。
KS検定をExcelを使って行う方法は下記の記事で紹介しました。
【Excelで簡単！】KS検定で出荷数が正規分布に従うかどうかを検定する
本記事では同じデータを使ってPythonでKS検定し、同じ結果になるかどうかを検証します。
そして、日々の出荷データが正規分布に従わないようなC商品でも、週でまとめれば正規分布になることをKS検定で示します。
PythonでKS検定する方法
まずはヒストグラムでばらつきをチェック
データは【Excelで簡単！】KS検定で出荷数が正規分布に従うかどうかを検定するで扱ったのと同じデータを使います。
Item Aが出荷数の多いA商品で、Item Bが出荷数の少ないC商品で、”Shipping_data_raw.csv”に入っています。
まずはヒストグラムでデータの分布をチェックしてみましょう。
Pythonで「ヒストグラム＋正規分布の近似曲線」を作成する４通りの方法でやった方法の中で、Pandasを使う方法でItem Aのヒストグラムを描いてみましょう。
import pandas as pd
import numpy as np
from scipy.stats import norm
data = pd.read_csv('Shipping_data_raw.csv')
xmin, xmax = data.Item_A.min(), data.Item_A.max()
mean, std = data.Item_A.mean(), data.Item_A.std()
x_axis = np.linspace(xmin, xmax, 100)
pdf = norm.pdf(x_axis, mean, std)
pdf_graph = pd.DataFrame({"X":x_axis, "norm":pdf}).set_index("X")
pdf_graph.plot()
data.Item_A.hist(bins=10, density=True)
同様にItem Bについては次のようになります。
Item Bは釣鐘状分布の右半分しかデータがなく、正規分布というには無理があることが一目でわかります。
Kstestを使って瞬時にKS検定
Item A
Scifyの中のstats.kstestモジュールを使えば1行でKS検定することができます。
kstest(data.Item_A, 'norm', args = (mean, std))kstestの後のカッコの中には
（検定したいデータのリスト、適合度を調べたい確率分布、確率分布のパラメータ）
を入力するだけです。
正規分布の場合のパラメータは平均と標準偏差です。
ヒストグラム＋正規分布近似曲線を描いた時に、既に平均と標準偏差をmeanとstdとして求めているのでそれを使いました。
実行結果は次のようになりました。
KstestResult(statistic=0.08357706947595453, pvalue=0.5213904505682965)KS値が0.0836、p値が0.521です。
KS値とは各データと正規分布との差の絶対値の最大値です。
これは【Excelで簡単！】KS検定で出荷数が正規分布に従うかどうかを検定するで既に求めていて、0.0849でした。
少し異なっていますね。
これは先ほど求めた標準偏差がサンプルの標準偏差であるためです。
Pandasで母集団の標準偏差である不偏標準偏差を求めるには、カッコの中で引数としてddof = 0を指定します。
std = data.Item_A.std(ddof = 0)このようにすることにより、stdに不変標準偏差が入ります。
まとめると、次のコードでKS検定できます。
mean, std = data.Item_A.mean(), data.Item_A.std(ddof = 0)
kstest(data.Item_A, 'norm', args = (mean, std))
KstestResult(statistic=0.08486080924359748, pvalue=0.5018970577950868)これでExcelで計算した場合のKS値と同じになりました。
またp値は0.5です。
これは0.05よりもはるかに大きい値なので、正規分布になるという帰無仮説を棄却できません。
正規分布になるといってよいでしょう。
ちなみに、正規分布のパラメータを求めるには、scify.statsのnorm.fit関数を使えばもっと簡単にできます。
para_A = norm.fit(data.Item_A)
kstest(data.Item_A, 'norm', args = para_A)これでも同じ結果になります。
Item B
次に同じようにしてItem BのKS検定も行ってみます。
para_B = norm.fit(data.Item_B)
kstest(data.Item_B, 'norm', args = para_B)
KstestResult(statistic=0.19572755079816284, pvalue=0.0015738427593380289)KS値はExcelで計算した結果と同じになりましたが、p値は0.0016と0.05より小さくなりました。
このことから、正規分布になるという帰無仮説は棄却されます。
出荷数の少ないItem Bの日々の出荷データは正規分布に従うとはいえません。
週単位の出荷数を計算する
しかし、日々の出荷数が正規分布にならないからといって、安全在庫理論が使えないというわけではありません。
週単位で出荷数をまとめ直してみましょう。
Item Bの出荷データは91日間のデータが並んだリストなので、これを頭から7つずつ取り出して足し合わせた新しいリストを作ってみましょう。
Pythonで行うには、リストの中の91個のデータをfor文で順に見ていきます。
その上で、0番目、7番目、14番目、、、というように7で割り切れる番数の時に、そこから6番先までのデータを足し合わせて、新しいリストに格納するという作業を繰り返せばよいでしょう。
例えば、次のようにコーディングできます。
Item_B = np.array(data['Item_B'])
Item_B_processed = []
cnt = 0
for i in Item_B:
if cnt % 7 == 0:
sum = np.sum(Item_B[cnt : cnt + 7])
Item_B_processed = np.append(Item_B_processed, sum)
cnt += 1
print(Item_B_processed)
[431. 541. 581. 425. 282. 459. 340. 510. 494. 639. 774. 330. 453.]これでItem Bの日単位の出荷データを週単位の出荷データに変換できました。
週単位の出荷データをヒストグラムにする
次に、週単位でまとめた出荷数がどのようなばらつき具合になっているかをヒストグラムで見てみましょう。
先ほどPandasのデータフレームをnumpy arrayに変換したので、そのデータをmatplotlibでヒストグラムにしてみます。
plt.hist(Item_B_processed, bins=10, density=True)
mean , std = Item_B_processed.mean(), Item_B_processed.std()
xmin, xmax = plt.xlim()
x_axis = np.linspace(xmin, xmax, 100)
pdf = norm.pdf(x_axis, mean, std)
plt.plot(x_axis, pdf)
plt.show()
だいぶ正規分布っぽくなってきました。
週単位の出荷データをKS検定する
それではこれが正規分布になっているといえるかどうかをKS検定してみましょう。
先ほどと全く同じようにやってみます。
para_B = norm.fit(Item_B_processed)
kstest(Item_B_processed, 'norm', args = para_B)
KstestResult(statistic=0.10785066751771366, pvalue=0.9939414472448639)するとp値は0.99となり、正規分布になるという帰無仮説を棄却できません。
正規分布に従うと仮定し、安全在庫理論を適用しても差し支えないでしょう。
まとめ
Pythonを使えば、KS検定が実質1行でできることがわかりました。
何千、何万というSKU数がある会社にとっては、Excelで分析するより遥かに短時間でデータ分析ができることになります。
また、日々の出荷データが少なすぎて正規分布にならない場合でも、週単位にまとめれば正規分布になる例を紹介しました。
実際には、この例よりも更に少ない出荷数の場合もあるでしょう。
しかしその場合でも、10日単位でまとめたり、1か月単位でまとめたりすることで正規分布になる可能性は大いにあるでしょう。
そのような解析にもPythonは大きな威力を発揮しそうです。
なお、ライブラリーを使えば今回のように１行でKS検定できてしまいますが、それでは何が行われているか原理がよくわからないので、こちらの記事でライブラリーを使わずにやってみました。
Pythonを使って経験分布と正規分布をグラフで視覚化してKS検定をしてみた