<jupyter_start><jupyter_text>顧客の評価 7.1 顧客のセグメンテーション 7.1.1 セグメンテーションの基準 7.1.2 クラスタ分析によるセグメンテーション<jupyter_code>import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from sklearn.linear_model import LinearRegression, LogisticRegression
from sklearn.metrics import confusion_matrix

%matplotlib inline
df_id_pos = pd.read_excel('978-4-274-22101-9.xlsx', 'ID付きPOSデータ（POSデータ）')
df_id_pos.head()
df_id_receipt = df_id_pos.groupby(['顧客ID', 'レシートNo'])['税抜価格'].sum()
loyal8 = df_id_receipt.count('顧客ID') >= 8
loyal8_cust_ids = loyal8[loyal8].index.values  # boolean array filtering

def tp(x):
    if x <=12:
        return '午前'
    elif x <= 15:
        return '午後'
    elif x <= 18:
        return '夕方'
    else:
        return '夜間'

df_id_pos0 = pd.DataFrame(df_id_pos)
df_id_pos0['時間帯'] = df_id_pos0['時間'].apply(tp)
df_id_pos0 = df_id_pos0[df_id_pos0['顧客ID'].isin(loyal8_cust_ids)]
df_id_tp_amount = pd.DataFrame(df_id_pos0.groupby(['顧客ID', '時間帯', 'レシートNo'])['税抜価格'].sum())
df_id_tp_amount.columns = ['購買金額']
df_id_tp_amount = df_id_tp_amount.mean(level=['顧客ID', '時間帯']).round()
df_id_tp_amount = df_id_tp_amount.unstack().fillna(0)
df_id_tp_amount
matplotlib.pyplot.rcParams['figure.figsize'] = (20.0, 5.0)

from scipy.cluster.hierarchy import linkage, dendrogram, fcluster

linkage = linkage(df_id_tp_amount, method='ward' ,metric='euclidean')
dendrogram(linkage)
plt.title("Dedrogram")
plt.ylabel("Threshold")
plt.show()

# t = 0.7 * max(linkage[:, 2])
t = 0.6 * max(linkage[:, 2])
group = fcluster(linkage, t, criterion='distance')
print(group)
df_id_tp_amount['group'] = group
# クラスタ1
df_id_tp_amount[df_id_tp_amount['group'] == 1].drop(['group'], axis=1).boxplot()
# クラスタ2
df_id_tp_amount[df_id_tp_amount['group'] == 2].drop(['group'], axis=1).boxplot()
# クラスタ3
df_id_tp_amount[df_id_tp_amount['group'] == 3].drop(['group'], axis=1).boxplot()
# クラスタ4
df_id_tp_amount[df_id_tp_amount['group'] == 4].drop(['group'], axis=1).boxplot()<jupyter_output><empty_output><jupyter_text>7.2 優良顧客の評価 7.2.1 RFM分析<jupyter_code>df_id_pos.head()
df_rfm = pd.DataFrame()

# １日〜25日までが集計対象、Rの基準日は26日
df_id_pos_25 = df_id_pos[df_id_pos['日'] <= 25]

df_rfm['R'] = 26 - df_id_pos_25.groupby(['顧客ID'])['日'].max()
df_rfm['F'] = df_id_pos_25.groupby(['顧客ID', '日']).size().groupby(['顧客ID']).size()
df_rfm['M'] = df_id_pos_25.groupby(['顧客ID'])['税抜価格'].sum()
df_rfm.head()
# ヒストグラムと散布図
pd.plotting.scatter_matrix(df_rfm)
# RFMそれぞれをクラス分けして購買生起確率を求める
cls = ['A', 'B', 'C']
keys= ['R' + r + '-F' + f + '-M' + m for r in cls for f in cls for m in cls]

conds = [
    [lambda c: c == 1, lambda c: (c >= 2) & (c <= 5), lambda c: c >= 6],
    [lambda c: c >= 7, lambda c: (c >= 3) & (c <= 6), lambda c: c <= 2],
    [lambda c: c >= 12551, lambda c: (c >= 4151) & (c <= 12550), lambda c: c <= 4150]
]

# df_rfm をクラス毎に分ける
df_rfm_classified = [df_rfm[conds[0][r](df_rfm['R']) & conds[1][f](df_rfm['F']) & conds[2][m](df_rfm['M'])] for r in range(3) for f in range(3) for m in range(3)]
# クラス毎の顧客IDをsetにする
rfm_classified_id_sets = map(lambda df: set(df.index.values), df_rfm_classified)
# 26日〜30日に購買があった顧客IDを抽出
target_id_set = set(df_id_pos[(df_id_pos['日'] >= 26) & (df_id_pos['日'] <= 30)]['顧客ID'].values)

props = []
for c in rfm_classified_id_sets:
    props.append(len(c & target_id_set) / len(c) if len(c) > 0 else 0)
plt.bar(keys, props)
plt.xticks(rotation='vertical')
plt.grid(True)
# R/F のクラスによる購買生起確率
# RとFだけ条件にして同様にクラス分け
df_rf_classified = [df_rfm[conds[0][r](df_rfm['R']) & conds[1][f](df_rfm['F'])] for r in range(3) for f in range(3)]
rf_classified_id_sets = map(lambda df: set(df.index.values), df_rf_classified)
props = []
for c in rf_classified_id_sets:
    props.append(len(c & target_id_set) / len(c) if len(c) > 0 else 0)

import numpy as np
pd.options.display.precision = 3

rf_keys = np.array([['R' + r, 'F' + f] for r in cls for f in cls])
df_rf_pp = pd.DataFrame({'prop': props})
df_rf_pp['R'] = rf_keys[:, 0]
df_rf_pp['F'] = rf_keys[:, 1]
df_rf_pp.set_index(['R', 'F']).unstack()<jupyter_output><empty_output><jupyter_text>7.2.2 ロジスティック回帰分析による来店予測<jupyter_code># 分析データ
df_rfm0 = df_rfm[df_rfm['M'] > 0]
df_rfm0['R2'] = df_rfm0['R'] ** 2
df_rfm0['M'] /= 1000
df_rfm0['再来店'] = list(map(lambda id: 1 if id in target_id_set else 0, df_rfm.index.values))
df_rfm0.head()
def std_t_p(X, Y, model):
    from scipy import stats
    
    params = np.append(model.intercept_, model.coef_)
    predictions = model.predict(X)
    newX = pd.DataFrame({ "Constant": np.ones(len(X)) }).join(pd.DataFrame(X))
    MSE = (sum((Y - predictions)**2)) / (len(newX) - len(newX.columns))

    var = MSE * (np.linalg.inv(np.dot(newX.T, newX)).diagonal())
    std = np.sqrt(var)
    t_value = params / std
    p_value =[2 * (1 - stats.t.cdf(np.abs(i), (len(newX) - 1))) for i in t_value]
    return std, t_value, p_value
X = df_rfm0.drop(['再来店'], axis=1).values
Y = df_rfm0['再来店'].values
model = LogisticRegression()
model.fit(X, Y)

std, t_value, p_value = std_t_p(X, Y, model)

df_lr_result = pd.DataFrame(index=['切片', 'R', 'F', 'M', 'R2'])
df_lr_result['偏回帰係数'] = np.append(model.intercept_, model.coef_)
df_lr_result['標準偏差'] = std
df_lr_result['t値'] = t_value
df_lr_result['p値'] = p_value
df_lr_result
# 予測来店確率
probabilities = model.predict_proba(X)[:, 1]
df_lr_prob = pd.DataFrame()
df_lr_prob['再来店なし'] = 1 - probabilities
df_lr_prob['再来店有り'] = probabilities
df_lr_prob.boxplot()
# 混同行列
cm = confusion_matrix(Y, predictions)
df_cm = pd.DataFrame(index=['実際来店有り', '実際来店なし'])
df_cm['予測来店有り'] = cm[:, 0]
df_cm['予測来店なし'] = cm[:, 1]
df_cm
# statsmodels で回帰
import statsmodels.api as st
X = df_rfm0.drop(['再来店'], axis=1).values
Y = df_rfm0['再来店'].values
X = st.add_constant(X, prepend = False)
mdl = st.Logit(Y, X)
mdl_fit = mdl.fit(method='lbfgs', maxiter=110)
mdl_fit.summary()<jupyter_output><empty_output><jupyter_text>7.3 因子分析・共分散構造分析による顧客の潜在的ニーズの構造分析<jupyter_code>df_survey_sheet = pd.read_excel('978-4-274-22101-9.xlsx', 'アンケート・データ')
df_survey_sheet
df_survey = df_survey_sheet.iloc[0:85, 1:11]
df_survey_x = df_survey.drop(['y1', 'y2'], axis=1)
df_survey_x
# Execute command below, restart kernel.
! conda install -y -c desilinguist factor_analyzer
from factor_analyzer import FactorAnalyzer

# 回転なし
model = FactorAnalyzer()
model.analyze(df_survey_x, 3, rotation=None)
model.loadings
model.get_factor_variance()
# varimax回転
model.analyze(df_survey_x, 3, rotation='varimax')
fa_varimax = model.loadings
fa_varimax
model.get_factor_variance()
# promax回転
model.analyze(df_survey_x, 3, rotation='promax')
fa_promax = model.loadings
fa_promax
model.get_factor_variance()<jupyter_output><empty_output><jupyter_text>7.3.2 共分散構造分析による因子間因果関係の分析<jupyter_code>df_survey_sheet = pd.read_excel('978-4-274-22101-9.xlsx', 'アンケート・データ')
df_survey = df_survey_sheet.iloc[0:85, 1:11]
df_survey.head()
# R 3.5 関連パッケージインストール
# sudo apt-key adv --keyserver keys.gnupg.net --recv-key 'E19F5F87128899B192B1A2C2AD5F960A256A04AF'
# sudo add-apt-repository 'deb http://cran.ism.ac.jp//bin/linux/debian stretch-cran35/'
# sudo apt update
# sudo apt install r-base, libxml2-dev
# mkdir ./rlib
# R -e "install.packages('XML', repos='http://cran.ism.ac.jp', lib=''./rlib')"
# R -e "install.packages('lavaan', repos='http://cran.ism.ac.jp', lib='./rlib')" >/tmp/R.install.packages.log1 2>&1
# R -e "install.packages(semPlot, repos=http://cran.ism.ac.jp, lib=./rlib)" >/tmp/R.install.packages.log 2>&1
! pip install pyper 
import pyper
df_survey_std = (df_survey - df_survey.mean()) / df_survey.std()
r = pyper.R(RCMD="/usr/bin/R", use_numpy='True', use_pandas='True')
r.assign("dat", df_survey_std)
r("source(file='./lavaan.R')")
r.get('rel')
print(r("summary(fit, standardized=TRUE, fit.measure=TRUE)"))<jupyter_output>try({summary(fit, standardized=TRUE, fit.measure=TRUE)})
lavaan 0.6-3 ended normally after 42 iterations

  Optimization method                           NLMINB
  Number of free parameters                         26

  Number of observations                            85

  Estimator                                         ML
  Model Fit Test Statistic                      33.445
  Degrees of freedom                                29
  P-value (Chi-square)                           0.260

Model test baseline model:

  Minimum Function Test Statistic              432.087
  Degrees of freedom                                45
  P-value                                        0.000

User model versus baseline model:

  Comparative Fit Index (CFI)                    0.989
  Tucker-Lewis Index (TLI)                       0.982

Loglikelihood and Information Criteria:

  Loglikelihood user model (H0)              -1001.747
  Loglikelihood unrestricted model (H1)       -985.025

  Number of fr[...]<jupyter_text>![上記モデルのパス図](./sem.png)  7.4 確率選択モデルによる購買行動モデル テキストで扱っているモデルはRのmlogitのドキュメントからすると conditional logit と呼ばれるもののようで、コレを計算するpythonモジュールはない模様。<jupyter_code>df_prod_select = pd.read_excel('978-4-274-22101-9.xlsx', '商品選択データ')
df_prod_select = df_prod_select.iloc[0:200, 0:9]
df_prod_select.head()
df_prod_select = df_prod_select.rename({'チラシD': 'チラシPB' }, axis=1)
df_prod_select.head()
# R版
# R -e "install.packages('mlogit', repos='http://cran.ism.ac.jp', lib='./rlib')"
# https://www.rdocumentation.org/packages/mlogit/versions/0.4-1/topics/mlogit

import pyper

# 日本語が使えない？ので名称をascii化する
# PBを基準にしたいのでA-Cの先頭にZをつけておく。アルファベット順？で先頭が基準になる模様
df_logit_e = pd.DataFrame()
df_logit_e['prod'] = df_prod_select['選択商品'].map({'A': 'ZA', 'B': 'ZB', 'C': 'ZC'})
df_logit_e['price.ZA'] = df_prod_select['価格A']
df_logit_e['price.ZB'] = df_prod_select['価格B']
df_logit_e['price.ZC'] = df_prod_select['価格C']
df_logit_e['price.PB'] = df_prod_select['価格PB']
df_logit_e['flyer.ZA'] = df_prod_select['チラシA']
df_logit_e['flyer.ZB'] = df_prod_select['チラシB']
df_logit_e['flyer.ZC'] = df_prod_select['チラシC']
df_logit_e['flyer.PB'] = df_prod_select['チラシPB']
r = pyper.R(RCMD="/usr/bin/R", use_numpy='True', use_pandas='True')
r(".libPaths('./rlib')")
r("library(mlogit)")
r.assign("df", df_logit_e)
r("data   <- mlogit.data(df, shape='wide', varying=2:9, choice='prod')")
head = r("head(data, 10)")
print(head)
r("result <- mlogit(prod ~ price + flyer, data)")
print(r("summary(result)"))<jupyter_output>try({summary(result)})

Call:
mlogit(formula = prod ~ price + flyer, data = data, method = "nr")

Frequencies of alternatives:
      PB       ZA       ZB       ZC 
0.000000 0.664557 0.297468 0.037975 

nr method
18 iterations, 0h:0m:0s 
g'(-H)^-1g = 9.46E-07 
gradient close to zero 

Coefficients :
                  Estimate  Std. Error z-value  Pr(>|z|)    
ZA:(intercept)  2.1009e+01  1.6955e+03  0.0124 0.9901137    
ZB:(intercept)  2.0153e+01  1.6955e+03  0.0119 0.9905168    
ZC:(intercept)  1.8776e+01  1.6955e+03  0.0111 0.9911644    
price          -3.3231e-02  5.7964e-03 -5.7329  9.87e-09 ***
flyer           1.8321e+00  4.9460e-01  3.7042 0.0002121 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Log-Likelihood: -89.323
McFadden R^2:  NaN 
Likelihood ratio test : chisq = NaN (p.value = NA)
{'fraction_non_alphanumeric': 0.1719304632747793, 'fraction_numerical': 0.04851187767361427, 'mean_word_length': 4.3731051344743275, 'pattern_counts': {'":': 1, '<': 35, '<?xml version=': 0, '>': 39, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 1}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '534284', 'license': 'MIT', 'n_tokens_mistral': 5458, 'n_tokens_neox': 5119, 'n_words': 832, 'path': 'notebooks/chapter7.ipynb', 'repo_name': 'ihiroky/statistics_for_marketing', 'size': '10987'}