アヒル本１０章多項ロジスティック回帰をpymcでやってみる。

8.4章のロジスティック回帰は、ある一つの項目（出席/欠席）についてのみ考えるが、多項ロジスティック回帰は複数の項目について考える。

StanとRでベイズ統計モデリング （Wonderful R） [ 松浦 健太郎 ]

カテゴリカル分布を使ってモデル化するのだが、各パラメータの１番目の要素を0に固定化するというのが味噌で、これをpymcで実装してみた。

# データの読み込み
df = pd.read_csv('data/data-category.txt')
# データの整形
X = df[['Age','Sex','Income']].values
Y = df.Y - 1
K = len(df.Y.unique()) # カテゴリ数

モデルの実装。

with pm.Model() as model:
    # パラメータ数は、後で0を加えるため、カテゴリ数-1になる。
    alpha = pm.Normal('alpha', mu=0, sd=100, shape=(K-1,))
    beta = pm.Normal('beta', mu=0, sd=100, shape=(3,K-1))

    # 各パラメータの1つを0に固定化する。
    alpha_f = pm.Deterministic('alpha_f', tt.concatenate([tt.zeros(1), alpha]))
    beta_f = pm.Deterministic('beta_f', tt.concatenate([tt.zeros((3,1)) , beta], axis=1))

    #
    mu = alpha_f + pm.math.dot(X, beta_f)
    theta = tt.nnet.softmax(mu)
    y = pm.Categorical('y', p=theta, observed=Y)

with model:
    start = pm.find_MAP(method='Powell')
    step = pm.NUTS()
    trace_sf = pm.sample(3000, step, start)

logp = -inf, ||grad|| = 0: 100%|██████████| 510/510 [00:00<00:00, 637.82it/s] 
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [beta, alpha]
Sampling 2 chains: 100%|██████████| 7000/7000 [03:05<00:00, 37.79draws/s]

MCMC実行結果。

pm.traceplot(trace_sf)

書籍には推定結果は載っていなかったため、データと推測結果を比較してみた。

# 推定結果の取り出し
alpha_f = np.mean(trace_sf['alpha_f'], axis=0)
beta_f = np.mean(trace_sf['beta_f'], axis=0)

各データで一番可能性が高いカテゴリを予測する。

preds = []
for uid in df.index:
    mu = alpha_f + np.dot(X[uid], beta_f)
    theta = tt.nnet.softmax(mu).eval()[0]
    argmax = np.argmax(theta)
    preds.append(argmax)

実際の結果と予測結果を比較する。0が最も多いはず。

(Y - preds).value_counts()

 0    101
 1     40
 2     32
-2     26
 4     23
-1     17
 3     16
 5     14
-4     13
-3     11
-5      7
Name: Y, dtype: int64

実際-予測=0が一番多かった。ということで、うまく推定できたと判断しよう。

以上