アイドルの人気投票における多項ロジスティック回帰の適用例

以前、アヒル本の多項ロジスティック回帰についてpymc3版を紹介した。

gaiasky.hatenablog.com

今回は、多項ロジスティック回帰の例として、「μ's とAqours の人気の差」を題材とした記事があったので、これを紹介したいと思う。これらの記事ではモデルはStanで実装されていたので、これをpymc3でトレースしてみることにする。

http://mikuhatsune.hatenadiary.com/entry/20170320/1490011326

http://abrahamcow.hatenablog.com/entry/2017/06/24/154902

データ生成のモデルは以下の通り

各9人が所属する2つのグループがある。
- メンバーの効果（人気）を $b_i$ 、μ'sグループの効果（人気）を $g$ とする。
- メンバー $i$ が所属するグループを $G_i$ とする。 $G_i=1$ はμ's、 $G_i=0$ はAqousに所属。
各メンバーの投票数とする。
- 得票数は多項分布に従う。 $v_i \sim Multinomial(p_i)$
- $p_i = softmax(\alpha_i)$
- $\alpha_i = g * G_i + b_i$

まずはデータの作成

dat_str = """name vote group
矢澤にこ 384 muse
西木野真姫 384 muse
南ことり 370 muse
東條希 336 muse
綾瀬絵里 283 muse
園田海未 263 muse
小泉花陽 251 muse
星空凛 226 muse
津島善子 190 aqours
黒澤ルビィ 171 aqours
国木田花丸 155 aqours
渡辺曜 141 aqours
高坂穂乃果 128 muse
松浦果南 110 aqours
黒澤ダイヤ 103 aqours
桜内梨子 93 aqours
小原鞠莉 80 aqours
高海千歌 74 aqours"""

df = pd.read_csv(io.StringIO(dat_str), sep=' ')

名前とグループをコード化する。

name_cate = pd.Categorical(df.name)
group_cate = pd.Categorical(df.group)
df['name_code'] = name_cate.codes
df['group_code'] = group_cate.codes
df.sort_values('name_code', inplace=True) # name_code順に並び替える

df.head()

	name	vote	group	name_code	group_code
2	南ことり	370	muse	0	1
10	国木田花丸	155	aqours	1	0
5	園田海未	263	muse	2	1
16	小原鞠莉	80	aqours	3	0
6	小泉花陽	251	muse	4	1

plt.figure(figsize=(9,6))
sns.barplot(x='vote', y='name', data=df, hue='group', orient='h')

f:id:gaiasky:20190908195638p:plain

pymc3によるモデルの実装とMCMC実行

with pm.Model() as model:
    # alpha[i] = g*G[i] + b[i] 
    # 各個人の効果
    b_mu = pm.Normal('b_mu', mu=0, sd=5)
    b_sd = pm.HalfCauchy('b_sd', beta=5)
    b = pm.Normal('b', mu=b_mu, sd=b_sd, shape=len(df)-1)
    b_f = pm.Deterministic('b_f', tt.concatenate([tt.zeros(1), b])) # 一人目の効果を0で固定
    # グループ効果
    g = pm.Normal('g', mu=0, sd=10)
    # 合計効果
    alpha = pm.Deterministic('alpha', g*df.group_code+b_f[df.name_code])
    # 
    p = pm.Deterministic('p', tt.nnet.softmax(alpha))
    y = pm.Multinomial('y', df.vote.sum(), p, observed=df.vote)

with model:
    trace = pm.sample(5000)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (2 chains in 2 jobs)
NUTS: [g, b, b_sd, b_mu]
Sampling 2 chains: 100%|██████████| 11000/11000 [01:05<00:00, 168.16draws/s]
The number of effective samples is smaller than 25% for some parameters.

pm.traceplot(trace, varnames=['g','b','p'])

f:id:gaiasky:20190908195834p:plain

推定結果の抽出

グループ（μ's）効果の確認。μ'sは、投票にに対してAquosより2.3倍の影響力がある。

pm.plot_posterior(np.exp(trace.get_values('g'))) #muse効果

f:id:gaiasky:20190908200030p:plain

可視化用のコード。

# 推定結果をグラフ表示するための整形関数
def mk_data4graph(df, trace, vname, do_exp):
    skip = 20
    if do_exp:
        estimate_params = np.exp(pm.trace_to_dataframe(trace, varnames=[vname]).loc[::skip]).T
    else:
        estimate_params = pm.trace_to_dataframe(trace, varnames=[vname]).loc[::skip].T        
    estimate_params.index = name_cate.categories
    
    # 平均値が高いキャラでソートする
    estimate_params['mean'] = estimate_params.mean(axis=1)
    estimate_params_sorted = estimate_params.sort_values('mean', ascending=False)
    estimate_params_sorted = estimate_params_sorted.drop(columns='mean')
    
    # グラフ用データに変換し、グループ名を追加
    data4graph = pd.merge(pd.melt(estimate_params_sorted.T), 
                          df[['name','group']], left_on='variable', right_on='name')
    return data4graph

各メンバーの投票率の可視化。

# 各メンバーの投票率）
plt.figure(figsize=(6,9))
sns.violinplot(x='value', y='variable', data=mk_data4graph(df, trace,'p', False),
               hue='group', palette=sns.color_palette("Set1", 2), orient='h',)
plt.xlabel('投票率$p_i$')
plt.title('メンバーの投票率$p_i$の推定値')

f:id:gaiasky:20190908200439p:plain

各メンバーの人気の可視化。

# 各メンバーの人気（exp(b)）
plt.figure(figsize=(6,9))
sns.violinplot(x='value', y='variable', data=mk_data4graph(df, trace,'b_f', True), 
               hue='group', palette=sns.color_palette("Set1", 2)[::-1], orient='h')
plt.xlabel('人気$\exp(b)$')
plt.title('メンバーの人気$\exp(b_i)$の推定値')