読者です 読者をやめる 読者になる 読者になる

tak0kadaの何でもノート

発声練習、生存確認用。

医学関連は 医学ノート

分布について - 「Think Stats」第4章

様々な分布の定義を確認して、実際のデータが分布に従っているか確認する方法を見ていく。連続分布に従って生成される実際の分布は離散的なので\(\rm PDF\)は使い勝手が悪い。\(\rm CDF\)について考える。

指数分布

\[ \begin{align} \rm (pd)f &= &\lambda e^{-\lambda x}\\ \rm (CD)F &= &1 - e^{-\lambda x}\\ mean(F) &= &1/\lambda\\ median(F) &= &\log(2)/\lambda \end{align} \]

である。

このCDFを変形してやると \[ log(CCDF) = log(1-CDF) = -\lambda x \] なのでCCDFをy軸をlogスケールにしてプロットすれば直線になるはずである。演習問題4-1に合わせて、平均32.6、44個の標本も含めてプロットすると

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_ylabel('ccdf')
ax.set_yscale('log')

for i, size in enumerate((44, 100, 1000, 100000)):
    pmf_5, bins = np.histogram(np.random.exponential(scale=32.6, size=size), bins=50, range=(0,250), density=True)
    ccdf = 1 - np.cumsum(pmf_5*5)
    bins = (bins + 2.5)[:-1]
    ax.plot(bins, ccdf, color=('k','r','g','b')[i], label=[str(j) for j in range(4)])

ax.legend([str(j)+'sample' for j in (44, 100, 1000, 100000)], loc='best')
plt.show()

f:id:kuyata:20140501184649p:plain

となる。

パレート分布

wikipediaによる式はかなりややこしい。簡略化されたものがテキストにはあって、 \[ \begin{align} \rm pdf &= &\dfrac{\alpha}{x_{m}}\biggl(\dfrac{x}{x_{m}}\biggr)^{-\alpha-1}\\ \rm cdf &= &1 - \biggl(\dfrac{x}{x_{m}}\biggr)^{-\alpha}\\ \end{align} \] ただし、\(x\)の最小値は\(x_{m}\)である。cdfを変形して、 \[ \log(\rm ccdf) = -\alpha(\log{x} - log{x_{m}}) \]

パレート分布は年収の分布を表すのに使えるらしい。京大卒の平均年収は677万円、先進国のパレート係数は1.5くらい(出典が分からない)というのが見えた。実際にシミュレーションしてみる。手元で計算してみたところ、\(x_{m}=400\)くらいになり、あまりにも現実にそぐわない気がするので撤退。

現実は以下

f:id:kuyata:20140501200716j:plain

参考: 全300校 出身大学別年収データ | 20代の”はたらき”データベース『キャリアコンパス』- powered by DODA -

パレート分布はスケールフリー性があるらしいのでプロットしてみる。

from __future__ import division
import numpy as np
from matplotlib import pylab

L = np.sort(np.random.pareto(1.5, size=100000))

pylab.scatter(L, L/np.sum(L))
pylab.show()

以下上から下へ順に拡大していっている。

f:id:kuyata:20140502001636p:plain

f:id:kuyata:20140502001646p:plain

f:id:kuyata:20140502001653p:plain

f:id:kuyata:20140502001701p:plain

正規分布

 \rm pdf = \frac{1}{\sqrt{2\pi\sigma^{2}}}exp\biggl(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\biggr)
 \rm cdf = \frac{1}{2}\biggl(1+erf\frac{x-\mu}{\sqrt{2\sigma^{2}}}\biggr)
\(\rm erf(x) = 2/\sqrt{\pi} \displaystyle\int_{0}^{x} e^{-t^{2}}\)

対数正規分布

\[ \rm cdf_{\rm lognormal} \equiv \rm cdf(log(x)) \]

対数正規分布の仕組みによると、むしろ日本の年収は対数正規分布に従うらしい。

from __future__ import division
import numpy as np
from matplotlib import pylab

pylab.hist(np.random.lognormal(mean=2,sigma=0.1,size=10000),bins=100, range=(0,20))
pylab.show()

salary = 300 + 10 * np.random.random_sample(1000)
pylab.hist(salary, range=(250,500), bins=500)
pylab.show()

for i in range(100):
    if i % 20 == 0:
        pylab.hist(salary, range=(250,1000), bins=500)
        pylab.show()

    up = []
    for j in 1000*np.random.random_sample(20):
        j = int(j)
        if j not in up:
            salary[j] = salary[j]*1.1
        up.append(j)

例に従ってプログラムを組んでみた結果は

f:id:kuyata:20140502004503p:plain

f:id:kuyata:20140502004516p:plain

f:id:kuyata:20140502004522p:plain

f:id:kuyata:20140502004530p:plain

f:id:kuyata:20140502004535p:plain

すごくギザギザな分布になった。どこかおかしいらしいがまたあす以降に検証したい。