<jupyter_start><jupyter_text>目的- エンゼルの出現確率を購入者毎に予測する- 独立に予測した結果と階層モデルで予測した結果を比較する<jupyter_code>import os, sys
import numpy as np
import pandas as pd
import scipy.stats as stats

import itertools
import math

%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns

import pymc3 as pm

sns.set()<jupyter_output><empty_output><jupyter_text>データの取得<jupyter_code>data = pd.read_csv('angel_data_for_publish.csv')
data.head()<jupyter_output><empty_output><jupyter_text>エンゼル出現確率の推定<jupyter_code>buyer = data['buyer'].values
buyer_idx = pd.Categorical(data['buyer']).codes
buyer_cat = pd.Categorical(data['buyer']).categories
data['buyer_idx'] = buyer_idx
lst_buyer = list(set(buyer_idx))
print(buyer_cat)
print(set(buyer_idx))
total_counts = data.groupby(['buyer_idx']).count()['silver'].values
angel_counts = data.query('silver > 0').groupby(['buyer_idx']).count()['silver'].values

print('total_count : {}'.format(total_counts))
print('angel_count : {}'.format(angel_counts))<jupyter_output>total_count : [142  80  59]
angel_count : [6 5 2]<jupyter_text>頻度（最尤推定）を計算<jupyter_code>theta_mle = angel_counts/total_counts
print(theta_mle)

fig = plt.figure(figsize=(8, 3))
ax = fig.subplots(1, 1)
cs = ['#FF4500', '#0000FF', '#00F1A1']
for idx in np.arange(0, len(theta_mle)):
    ax.vlines(theta_mle[idx], 0, 1, colors=cs[idx%len(cs)], label=buyer_cat[idx])
ax.set_xlim((0.0, 0.1))
ax.set_xlabel('$\\theta_{MLE}$')
ax.legend()
plt.savefig('buyer_effect_mle.png')<jupyter_output>[0.04225352 0.0625     0.03389831]<jupyter_text>購入者ごとに独立に推定<jupyter_code>with pm.Model() as comparing_buyer_i:
    theta = pm.Beta('theta', alpha=1, beta=1, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    
    trace_iso = pm.sample(5000, chains=1, random_seed=100)
pm.traceplot(trace_iso)
pm.summary(trace_iso)
fig = plt.figure(figsize=(8, 4))
ax = fig.subplots(1, 1)

for i in np.arange(len(buyer_cat)):
    sns.distplot(trace_iso['theta'][:,i], label=buyer_cat[i], ax=ax)
ax.legend()
ax.set_xlabel('angel rate')
ax.set_ylabel('frequent')
plt.savefig('buyer_effect_iso.png')
# 差の分布を確認
n_c = len(list(itertools.combinations(np.arange(0, len(lst_buyer)), 2)))
n_col = 3
n_row = math.ceil(n_c/n_col)
fig = plt.figure(figsize=(12, 3*n_row))
ax = fig.subplots(n_row, n_col)
if n_row == 1:
    ax = ax[np.newaxis, :]

cnt=0
for (i,j) in itertools.combinations(np.arange(0, len(lst_buyer)), 2):
    theta_diff = trace_iso['theta'][:, i] - trace_iso['theta'][:, j]
    pm.plot_posterior(theta_diff, ref_val=0, ax=ax[int(cnt/n_col), int(cnt%n_col)])
    ax[int(cnt/n_col), int(cnt%n_col)].set_title('{}-{}'.format(buyer_cat[i], buyer_cat[j]))
    cnt+=1
plt.savefig('buyer_effect_diff_isomodel.png')<jupyter_output><empty_output><jupyter_text>購入者の性質を加える- 真の出現率は決まっているはず(真の確率を$p$とする)- 購入者毎の運の要素が入ってくる（不正行為かも、独自の購入戦略があるのかも）可能性がある  - 購入者毎の特性を個人差$u_i$とする- エンゼルの出現は確率$\theta_i$の二項分布- 確率$\theta_i$は以下のロジットリンク関数で線形モデルを仮定$$logit(\theta_i) = p + u_i$$- $p$、$u_i$は正規分布を事前分布とする<jupyter_code>with pm.Model() as comparing_buyer_m1:
    su = pm.HalfNormal('su', sd=10)
    p = pm.Normal('p', mu=0, sd=10)
    u = pm.Normal('u', mu=0, sd=su, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', 
                        n=total_counts[lst_buyer], 
                        p=pm.math.sigmoid(p+u[lst_buyer]), observed=angel_counts[lst_buyer])
    
    trace_h1 = pm.sample(3000, chains=3, random_seed=100)
pm.traceplot(trace_h1)
pm.summary(trace_h1)
pm.plot_posterior(trace_h1, varnames=['p', 'u'])
fig = plt.figure(figsize=(8, 4))
ax = fig.subplots(1, 1)

for i in np.arange(len(buyer_cat)):
    sns.distplot(trace_h1['u'][:,i], label=buyer_cat[i], ax=ax)
ax.legend()
ax.set_xlabel('angel rate')
ax.set_ylabel('frequent')
plt.savefig('buyer_effect_h1model.png')
# 差の分布を確認
n_c = len(list(itertools.combinations(np.arange(0, 3), 2)))
n_col = 3
n_row = math.ceil(n_c/n_col)
fig = plt.figure(figsize=(12, 3*n_row))
ax = fig.subplots(n_row, n_col)
if n_row == 1:
    ax = ax[np.newaxis, :]

cnt=0
for (i,j) in itertools.combinations(np.arange(0, 3), 2):
    u_diff = trace_h1['u'][:, i] - trace_h1['u'][:, j]
    pm.plot_posterior(u_diff, ref_val=0, ax=ax[int(cnt/n_col), int(cnt%n_col)])
    ax[int(cnt/n_col), int(cnt%n_col)].set_title('{}-{}'.format(buyer_cat[i], buyer_cat[j]))
    cnt+=1
plt.savefig('buyer_effect_diff_h1model.png')
def sigmoid(x):
    return 1.0 / (1.0 + np.exp(-x))
post_min = pm.summary(trace_h1).loc['p', 'hpd_2.5']
post_mean = pm.summary(trace_h1).loc['p', 'mean']
post_max = pm.summary(trace_h1).loc['p', 'hpd_97.5']
print('{} < {} < {}'.format(sigmoid(post_min), sigmoid(post_mean), sigmoid(post_max)))
fig = plt.figure(figsize=(8, 4))
ax = fig.subplots(1, 1)

sns.distplot(trace_h1['p'], ax=ax)
ax.set_xlabel('angel rate')
ax.set_ylabel('frequent')
plt.savefig('buyer_effect_h1model_tp.png')<jupyter_output>No handles with labels found to put in legend.<jupyter_text>階層モデルとして、全体を統括するパラメータがあると仮定- 真の出現率を決めるパラメータがあるはず- 購入者毎の運によって↑のパラメータに基づいたパラメータがサンプルされる- エンゼルの出現は確率$\theta_i$の二項分布- $\theta$の事前分布はベータ分布- ベータ分布のパラメータに事前分布を設定<jupyter_code>with pm.Model() as comparing_buyer_m2:
    alpha = pm.HalfCauchy('alpha', beta=10)
    beta = pm.HalfCauchy('beta', beta=10)
    
    theta = pm.Beta('theta', alpha=alpha, beta=beta, shape=len(set(buyer_idx)))
    
    angel = pm.Binomial('angel', n=total_counts[lst_buyer], p=theta[lst_buyer], observed=angel_counts[lst_buyer])
    
    trace_h2 = pm.sample(2000, chains=3, random_seed=100)
pm.traceplot(trace_h2)
fig = plt.figure(figsize=(8, 4))
ax = fig.subplots(1, 1)

for i in np.arange(len(buyer_cat)):
    sns.distplot(trace_h2['theta'][:,i], label=buyer_cat[i], ax=ax)
ax.legend()
ax.set_xlabel('angel rate')
ax.set_ylabel('frequent')
plt.savefig('buyer_effect_h2model.png')
# 差の分布を確認
n_c = len(list(itertools.combinations(np.arange(0, len(lst_buyer)), 2)))
n_col = 3
n_row = math.ceil(n_c/n_col)
fig = plt.figure(figsize=(12, 3*n_row))
ax = fig.subplots(n_row, n_col)
if n_row == 1:
    ax = ax[np.newaxis, :]

cnt=0
for (i,j) in itertools.combinations(np.arange(0, len(lst_buyer)), 2):
    theta_diff = trace_h2['theta'][:, i] - trace_h2['theta'][:, j]
    pm.plot_posterior(theta_diff, ref_val=0, ax=ax[int(cnt/n_col), int(cnt%n_col)])
    ax[int(cnt/n_col), int(cnt%n_col)].set_title('{}-{}'.format(buyer_cat[i], buyer_cat[j]))
    cnt+=1
plt.savefig('buyer_effect_diff_h2model.png')
# thetaの事前分布
x = np.linspace(0, 1, 100)
for i in np.random.randint(0, len(trace_h2), size=100):
    pdf = stats.beta(trace_h2['alpha'][i], trace_h2['beta'][i]).pdf(x)
    plt.plot(x, pdf,  'C1', alpha=0.2)

dist = stats.beta(trace_h2['alpha'].mean(), trace_h2['beta'].mean())
pdf = dist.pdf(x)
mode = x[np.argmax(pdf)]
mean = dist.moment(1)
plt.plot(x, pdf, label='mode = {:.2f}\nmean = {:.2f}'.format(mode, mean))

plt.legend(fontsize=14)
plt.xlabel('$\\theta_{prior}$', fontsize=16)
plt.tight_layout()
plt.savefig('buyer_effect_prior_h2model.png', dpi=300, figsize=(5.5, 5.5))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1871004404034664, 'fraction_numerical': 0.030117914476488136, 'mean_word_length': 6.154471544715447, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '74598', 'license': 'MIT', 'n_tokens_mistral': 3545, 'n_tokens_neox': 3393, 'n_words': 436, 'path': 'analysis/jupyter/estimate-angel-rate-multi-pub.ipynb', 'repo_name': 'tok41/choco-ball-statistics', 'size': '7039'}