ポワソン分布のシミュレーション
ポワソン分布は単位時間あたりλ回起こる出来事の確率分布と説明されることが多い。
これをPythonのシミュレーションで確かめてみる。
In [373]:
import random
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pandas import DataFrame, Series
from scipy.stats import poisson
%matplotlib inline
シミュレーションを行うには単位時間あたりの長さのほかにシミュレートする時間の長さが必要だ。
単位時間あたり4回起こる出来事を長さが2の時間の中にランダムに配置するとすると、8個のイベントを配置するということになる。
In [374]:
[random.random() * 2 for _ in range(0, 8)]
Out[374]:
単位時間1の中にいくつ入ったのかを数える。
In [375]:
[random.random() * 2 < 1 for _ in range(0, 8)].count(True)
Out[375]:
これを10000回繰り返してヒストグラムを描く。
In [376]:
simulations = [
[random.random() * 2 < 1 for _ in range(0, 8)].count(True)
for _ in range (0, 10000)]
plt.hist(simulations, bins=range(0, 8+2), normed=True, color='gray')
plt.show()
これは実質的には試行回数が8で確率が1/2の二項分布のシミュレーションを行っていることになる。 8個のイベントしか配置していないから8回よりも多くなることはない。
一方ポワソン分布のパラメタは単位時間あたりの期待値だけで、8のような上限はない。 そこで期待値が4に固定されるようにして、長さとイベント数を増やしてみる。
さらにシミュレーション結果の隣にポワソン分布の理論値のヒストグラムを並べるようにする。
In [377]:
lam = 4
lengths = [2, 4, 8, 16]
fig_columns = 2
fig_rows = len(lengths) / fig_columns + (1 if len(lengths) % fig_columns > 0 else 0)
fig, axes = plt.subplots(fig_rows, fig_columns)
for i, length in enumerate(lengths):
num_events = int(length * lam) # 単位長さ当たりのイベント数num_events/length=lamになるようにする
simulations = [
# 長さlengthの中にnum_events個のイベントを配置し、1未満に入ったイベントの数を数える
[random.random() * length < 1 for _ in range(0, num_events)].count(True)
for _ in range (0, 10000)] # 10000回シミュレートする
ax = axes[i / fig_columns, i % fig_columns]
x_range = Series(range(0, int(poisson.ppf(0.99999, lam))))
(_, _, bar_sim) = ax.hist(simulations, bins=x_range, normed=True, color='gray', width=0.5)
bar_poi = ax.bar(x_range + 0.5, poisson.pmf(x_range, lam), color='white', width=0.5)
ax.set_xlabel('count')
ax.set_ylabel('probability')
ax.set_title('length=%d, # of events=%d' % (length, num_events))
fig.legend((bar_sim[0], bar_poi[0]), ('simulation', 'poisson'), 'upper right', bbox_to_anchor=(1.2, 1))
fig.subplots_adjust(wspace=0.5, hspace=1)
plt.show()
長さとイベント数が多くなるにつれてシミュレーションの分布がポワソン分布に近づいていることがわかる。
参考にした文献：「ポワソン分布・ポワソン回帰・ポワソン過程」島谷健一郎（近代科学社）
複数の観察データからのパラメタ推定 (PyMC3, PyStan)
今度は二項分布のパラメタが1つだけあって、同じ二項分布から複数のセッション（n回中k回）がデータとして取られる場合に、 その1つだけのパラメタを推定するということを考える。
In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
%matplotlib inline
PyMC3¶
In [2]:
import pymc3 as pm
from pymc3 import Beta, Binomial, Bernoulli, Model, Deterministic
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
PyMC3でモデルを定義する。 観察されるデータはk1をk2の2つあり、同じパラメタthetaの二項分布から生み出されているとする。2つの観察データは10回中5回と10回中7回であった。
In [3]:
model = Model()
with model:
# Prior
theta = Beta('theta', alpha=1, beta=1)
# Observed
k1 = Binomial("k1", n=10, p=theta, observed=5)
k2 = Binomial("k2", n=10, p=theta, observed=7)
trace = sample(5000)
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
上記ではk1とk2という2つのBinomialを記述したが、下記のように1つのBinomialでobservedを配列にする書き方もできるようだ。
In [4]:
model = Model()
with model:
# Prior
theta = Beta('theta', alpha=1, beta=1)
# Observed
k = Binomial("k", n=10, p=theta, observed=[5, 7])
trace = sample(5000)
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
ところで10回中5回、7回の2つの二項分布のセッションがあるというのは、20個のベルヌーイ分布で12回が1であるという見方をしてもいいのだろうか。
In [5]:
model = Model()
with model:
# Prior
theta = Beta('theta', alpha=1, beta=1)
# Observed
k = Bernoulli("k", p=theta, observed=[1,1,1,1,1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0])
trace = sample(5000)
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
この場合でもパラメタの推定値が同様であることがわかる。
PyStan¶
同じことをPyStanでやってみる。まずはk1とk2を別々に定義する書き方から。
In [6]:
import pystan
In [7]:
code = """
data {
int<lower=0> k1;
int<lower=1> n1;
int<lower=0> k2;
int<lower=1> n2;
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1, 1);
k1 ~ binomial(n1, theta);
k2 ~ binomial(n2, theta);
}
"""
data = {'k1': 5, 'n1': 10, 'k2': 7, 'n2': 10}
In [8]:
fit = pystan.stan(model_code=code, data=data)
print(fit)
Stanでも配列を使うことができる。forを使ってインデクシングするようだ。
In [9]:
code = """
data {
int N;
int<lower=0> k[N];
int<lower=1> n[N];
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1, 1);
for (i in 1:N)
k[i] ~ binomial(n[i], theta);
}
"""
data = {'N': 2, 'k': [5, 7], 'n': [10, 10]}
fit = pystan.stan(model_code=code, data=data)
print(fit)
2つの二項分布のパラメタの差の推定(PyMC3, PyStan)
二項分布から生成されると仮定される2つのデータがあったときにその差（2つの二項分布のパラメタの差）を推定することを考える。
In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
%matplotlib inline
PyMC3¶
In [2]:
import pymc3 as pm
from pymc3 import Beta, Binomial, Model, Deterministic
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
PyMC3でモデルを定義する。ここでは両方とも10回の施行で、それぞれ5回と7回が観察されたものとする。 パラメタの差を表現するにはDeterministicを使う。
In [3]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=10, p=theta1, observed=5)
k2 = Binomial("k2", n=10, p=theta2, observed=7)
MCMCする。
In [4]:
with model:
trace = sample(5000)
結果を表示する。
In [5]:
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
同じことを違うデータでやってみる。試行回数と観察回数をそれぞれ10倍する。
In [6]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=100, p=theta1, observed=80)
k2 = Binomial("k2", n=100, p=theta2, observed=70)
trace = sample(5000)
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
事後分布の左右の幅が狭くなったのがわかる。
次に試行回数が異なってそれぞれゼロ回の観察だった場合。
In [7]:
model = Model()
with model:
# Prior
theta1 = Beta('theta1', alpha=1, beta=1)
theta2 = Beta('theta2', alpha=1, beta=1)
# Difference
delta = Deterministic('delta', theta1 - theta2)
# Observed
k1 = Binomial("k1", n=1, p=theta1, observed=0)
k2 = Binomial("k2", n=5, p=theta2, observed=0)
trace = sample(5000)
summary(trace[1000:], varnames=['delta'])
traceplot(trace[1000:], varnames=['delta'])
plt.show()
両方とも0なので、差が0のところが最もありそうなのだが、k1の試行回数が少ない分確信がないので右の裾が長い。
PyStan¶
同じことをPyStanでやってみる。
In [8]:
import pystan
In [9]:
code = """
data {
int<lower=0> k1;
int<lower=1> n1;
int<lower=0> k2;
int<lower=1> n2;
}
parameters {
real<lower=0,upper=1> theta1;
real<lower=0,upper=1> theta2;
}
transformed parameters {
real<lower=-1,upper=1> delta;
delta <- theta1 - theta2;
}
model {
theta1 ~ beta(1, 1);
theta2 ~ beta(1, 1);
k1 ~ binomial(n1, theta1);
k2 ~ binomial(n2, theta2);
}
"""
data = {'k1': 5, 'n1': 10, 'k2': 7, 'n2': 10}
In [10]:
fit = pystan.stan(model_code=code, data=data)
In [11]:
print(fit)
fit.plot()
Out[11]:
PyMC3とPyStanによる二項分布のパラメタ推定
前回の記事で二項分布のパラメタ推定をPyMC2で行った。 しかしPyMC2をいまから使っていくのも微妙な気がしてきたので新しいPyMC3で書き直す。 ついでにPyStanでも同じことをして比べてみる。
In [1]:
import numpy as np
import pandas as pd
from pandas import DataFrame, Series
from matplotlib import pyplot as plt
%matplotlib inline
PyMC3¶
まずはPyMC3のモジュールをインポートする。 theano.config.warn.round=Falseは警告を抑止するための設定。
In [2]:
import pymc3 as pm
from pymc3 import Beta, Binomial, Model
from pymc3 import traceplot, sample, summary
import theano
theano.config.warn.round=False
次にモデルを定義する。
In [3]:
model = Model()
with model:
theta = Beta('theta', alpha=1, beta=1)
k = Binomial("k", n=10, p=theta, observed=5)
PyMC3ではこのようにwith構文を使ってモデルを定義する。 何をしているのかわからなくて面食らうが、見た目上モデル定義がまとまりをなすようにするための工夫だろうか。 多分ダイナミックな文脈を作ってBetaやBinomialなどが現在の文脈と結びつくようにしているのかもしれない。
MCMCを行うにはsample関数を使う。
In [4]:
with model:
trace = sample(5000)
summary関数やtraceplot関数で結果を表示することができる。
In [5]:
summary(trace[1000:])
traceplot(trace[1000:])
plt.show()
PyStan¶
次に同じことをPyStanでやってみる。
In [6]:
import pystan
PyStanではモデルは完全にStanの文法で記述することになる。
In [7]:
code = """
data {
int<lower=0> k;
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1, 1);
k ~ binomial(10, theta);
}
"""
data = {'k': 5}
stan関数でMCMCが行われる。
In [8]:
fit = pystan.stan(model_code=code, data=data)
stan関数の戻り値をprintしたりplotメソッドを呼ぶと結果を表示できる。
In [9]:
print(fit)
fit.plot()
Out[9]:
なんで2つずつ出るのだろうか。これはよくわからない。
PyMCの初歩
PyMCを試すことにする。 Pythonで体験するベイズ推論 を買った関係上、PyMC3ではなくてPyMC2を使っている。
この記事の内容は書籍と関係なく、書籍よりずっと初歩的なものである。
In [1]:
import pymc as pm
from pymc.Matplot import plot as mcplot
from pymc import Beta, Binomial, Uniform, Model, MCMC
from pymc.Matplot import plot as mcplot
from matplotlib import pyplot as plt
%matplotlib inline
データが二項分布によって生成されるモデルを考える。 二項分布では事象発生の確率を示すパラメタ$\theta$があって、これが推定の対象となる。
$\theta$は0から1までの連続値であり、特に事前の情報を持たないので、 0から1までのすべてが等確率（密度）となる事前分布を使う。 これは具体的にはパラメタ$\alpha$と$\beta$がそれぞれ1のベータ分布である。
学習に使うデータとして、10回中5回事象が起こったとしよう。
In [2]:
theta = Beta('$\\theta$', alpha=1, beta=1)
k = Binomial("k", n=10, p=theta, value=5, observed=True)
model = Model([theta, k])
このモデルを使ってMCMCを行い、結果をプロットする。
当然だが右側のプロットで示されるように、推定値は0.5の周りになる。 95%区間は点線で示されている。
10回中5回でなく、100回中50回だったら0.5に近いことにもっと確証を持つだろう。
In [4]:
95%区間が狭くなったことがわかる。
1回中1回、1回中0回、2回中1回も試してみる。
In [5]:
In [6]:
In [7]:
ところで二項分布のパラメタの事前分布としてベータ分布を使うというのは殆どの文献で慣例的にそうなっているのだが、 等確率分布ならばUniformを使ってもよいのだろうか。
In [8]:
In [9]:
結果に特に違いはないように見える。
In [ ]:
「データ指向のソフトウェア品質マネジメント」の「テスト工程での不具合検出数の予測」をPythonで行う
In [1]:
import math
import pandas as pd
from pandas import DataFrame, Series
import matplotlib.pyplot as plt
%matplotlib inline
まずは対象データを読み込む。これは書籍のサポートページからダウンロードできる。 ダウンロードに書籍記載のパスワードが必要なので生データをここには載せないが、 データ自体はFenton (2008)のものである。
ソースコード行数KLOCと不具合検出数Failures以外はカテゴリーデータ（順序尺度）で、D1は開発メンバーの経験、P5はステークホルダーの関与の度合い、P6は顧客とのコミュニケーション、P7は協力会社のマネジメントを表す。
In [2]:
data1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj')
data1.columns
Out[2]:
In [3]:
def scatter_corr_matrix(df):
pd.scatter_matrix(df)
plt.show()
return df.corr()
scatter_corr_matrix(data1)
Out[3]:
In [4]:
import statsmodels.formula.api as smf
smf.ols(formula="Failures ~ KLOC", data=data1).fit().summary()
Out[4]:
In [5]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC)", data=data1).fit().summary()
Out[5]:
In [6]:
data1_1 = pd.read_excel(u"4.2節/4-2-data.xls", index_col='Prj', sheetname='data1-1')
data1_1.describe()
Out[6]:
ソースコード行数がすでに説明変数に含まれているため、 不具合検出密度と各順序尺度の相関（スピアマンの順位相関係数）を計算する。
In [7]:
from scipy.stats import spearmanr
data1_1_d = data1_1.assign(Density = data1_1['Failures'] / data1_1['KLOC'])
data1_1_d.corr(method='spearman')
Out[7]:
不具合検出密度と各順序尺度の関係を箱ひげ図で表す。
In [8]:
data1_1_d.boxplot(column='Density', by='P7', return_type='dict')
data1_1_d.boxplot(column='Density', by='D1', return_type='dict')
data1_1_d.boxplot(column='Density', by='P5', return_type='dict')
data1_1_d.boxplot(column='Density', by='P6', return_type='dict')
plt.show()
これを見ると不具合検出密度と相関が高いのはP7（協力会社のマネジメントが高いと不具合検出密度が低い）だが、「データ志向の～」ではP7はデータ件数が少ない(n=7)こと、ソースコード行数との相関係数も極めて大きいこと（多重共線性の疑い）を理由に退けている。
重回帰分析の説明変数をステップワイズで選択¶
P7以外の順序尺度をステップワイズで選択する。まずは全部を使う。
In [9]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6", data=data1).fit().summary()
Out[9]:
ここから1つずつ減らしていく。
In [10]:
formulas = [
"np.log(Failures) ~ np.log(KLOC) + D1 + P5 + P6",
"np.log(Failures) ~ np.log(KLOC) + D1 + P5",
"np.log(Failures) ~ np.log(KLOC) + D1",
"np.log(Failures) ~ np.log(KLOC)",
]
steps = [[res.llf, res.aic, res.rsquared, res.rsquared_adj, formula]
for (res, formula) in [(smf.ols(formula=formula, data=data1).fit(), formula) for formula in formulas]]
DataFrame(steps, columns=['Log-Likelihood', 'AIC', 'R-squared', 'Adj. R-squared', 'formula'])
Out[10]:
AICおよびAdj. R-squaredで判断するとlog(KLOC)とD1を使うモデルが最も良い。
In [11]:
smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1", data=data1).fit().summary()
Out[11]:
このサマリーを見るとD1のダミー変数の中でも2_Lと3_MのP値が高い。そこで「データ指向の～」では3_M以下の水準を統合することを提案している。
In [12]:
d1_map = {'1_VL': '3_M', '2_L': '3_M', '3_M': '3_M', '4_H': '4_H', '5_VH': '5_VH'}
data1_d1_new = data1.assign(D1_new = data1['D1'].apply(lambda x: d1_map[x]))
ols_final = smf.ols(formula="np.log(Failures) ~ np.log(KLOC) + D1_new", data=data1_d1_new).fit()
ols_final.summary()
Out[12]:
In [13]:
scatter_corr_matrix(
pd.concat(
[data1.dropna(subset=['KLOC']).reset_index(drop=True),
Series(np.exp(ols_final.predict()), name='predict')],
axis=1)[['Failures', 'predict']])
Out[13]:
上記では散布図も表示している。なお、「データ指向の～」（2015年1月13日第4刷）に掲載の散布図はX軸とY軸が逆のように思える（軸ラベル上Xが予測でYが実績なのだが、右上でなく左下の散布図になっている）。
相関係数が0.887であることから、うまく予測できているように見えると判断している。
このくだりにはちょっとよくわからないことろがあって、そもそも予測値と実績値の相関係数というのは予測モデルの評価の仕方として適当なのだろうかという疑問がある。 というのはモデルが誤っていて「本来の値の約半分の値を予測する」というものであった場合、予測モデルとしては役に立たないが、相関係数は半分でない場合と同様に高くなるはずだからだ。
次に予測誤差 MRE (Megnitude of Relative Error)と呼ばれる評価指標を出している。こちらのほうが直感的な定義に思える。MREは1つ1つの予測に対する値なので、全体で1つの指標値にするには予測誤差の平均とか中央値とかで見るようだ。
In [14]:
def mre(actual, predicted):
df = DataFrame({'actual': actual, 'predicted': predicted})
mre = df.eval("abs(actual - predicted) / actual")
mre.plot.box()
plt.show()
print mre.describe().round(2)
return mre.median()
mre(
data1.dropna(subset=['KLOC'])['Failures'],
np.exp(ols_final.predict()))
Out[14]:
In [16]:
import statsmodels.api as sm
poisson_model = smf.glm(
formula='Failures ~ np.log(KLOC) + D1 + P5 + P6',
data=data1,
family=sm.families.Poisson()).fit()
poisson_model.summary()
Out[16]:
In [17]:
scatter_corr_matrix(
pd.concat(
[data1.dropna(subset=['KLOC']).reset_index(drop=True),
Series(poisson_model.predict(), name='predict')],
axis=1)[['Failures', 'predict']])
Out[17]:
予測値と実測値の相関係数は向上している。
一方MREを計算すると中央値では改善するが平均値では悪くなっており、またばらつきも大きくなっている。
In [18]:
mre(
data1.dropna(subset=['KLOC'])['Failures'],
poisson_model.predict())
Out[18]:
In [ ]:
カテゴリーデータを説明変数に含む回帰分析
今度は説明変数にカテゴリーデータを含む場合の回帰分析を行う。
In [77]:
%matplotlib inline
対象となるデータはR処理系のcarパッケージに付属しているPrestigeというデータを write.csv(Prestige, "Prestige.csv", quote=FALSE, row.names=TRUE) でCSVにしたものである。
このデータの各行は職業で、カナダの国勢調査で得られたデータである。変数は「定性的データ分析」（金明哲）によると次の通り。
education: 平均教育年数, income: 平均所得, women: 女性が占める割合, prestige: 職業地位の得点, census: カナダの国勢調査の職業コード, type: {bc: ブルーカラー, prof: 専門家, wc: ホワイトカラー}
職業コードは無視すると、typeがカテゴリーデータである。
次に散布図行列を書く。pandasのscatter_matrixは残念ながらカテゴリーデータを散布図行列に含めてくれない。もっとも含まれていても相関の解釈はしづらい。
In [79]:
from pandas.tools.plotting import scatter_matrix
scatter_matrix(prestige.drop("census", axis=1), figsize=(7, 7))
plt.show()
目的変数にしたいprestigeと相関を持ちそうなのがeducationとincomeであることがわかる。
重回帰分析にはstatsmodels.formula.apiのほうの関数を使う。というのもこちらを使うとカテゴリー変数に対して自動的にダミー変数を作ってくれるからだ。
In [80]:
import statsmodels.formula.api as smf
mod = smf.ols(formula="prestige ~ education + income + women + type", data=prestige)
res = mod.fit()
print res.summary()
次に変数選択を行う。ここではAICでもAdj. R-squaredでも同じく、education + income + typeの場合が最もよい。
In [81]:
formulas = [
"prestige ~ education + income",
"prestige ~ education + income + type",
"prestige ~ education + income + women + type",
]
steps = [[res.llf, res.aic, res.rsquared, res.rsquared_adj, formula]
for (res, formula) in [(smf.ols(formula=formula, data=prestige).fit(), formula) for formula in formulas]]
DataFrame(steps, columns=['Log-Likelihood', 'AIC', 'R-squared', 'Adj. R-squared', 'formula'])
Out[81]:
選択された変数で改めて重回帰分析を行う。
In [82]:
mod = smf.ols(formula="prestige ~ education + income + type", data=prestige)
res = mod.fit()
print res.summary()
多重共線性について確かめるためVIFを計算する。
In [83]:
from statsmodels.stats.outliers_influence import variance_inflation_factor
num_cols = mod.exog.shape[1] # 説明変数の列数
vifs = [variance_inflation_factor(mod.exog, i)
for i in range(0, num_cols)]
DataFrame(vifs, index=mod.exog_names, columns=["VIF"])
Out[83]:
ここではカテゴリー変数typeのダミー変数に対してVIFが計算されているが、 「定性的データ分析」のRの実行結果ではtypeそのものに対して6.102131という結果が出ている。 これをStatsmodelsでどうすればいいのか（Rがどう計算しているのか）はわからない。 いったん10以下であることから多重共線性はないと判断する。
「定性的データ分析」ではprestigeに対するincomeの関係が直線でなく右側が吊り上がっていることに注目し、変数incomeを対数変換することを提案している。 これをStatsmodelsで行う場合、式にはlog(income)でなく、numpy.logを使って書く必要があるようだ。
In [84]:
import math
mod = smf.ols(formula="prestige ~ education + np.log(income) + type", data=prestige)
res = mod.fit()
print res.summary()
AICが良くなり、Adj. R-squaredも0.828から0.849に向上している。
