A群（Control）[465.0, 458.6, 466.5, 475.2, 457.7, 457.7, 475.8, 467.7, 455.3, 465.4, 455.4, 455.3, 462.4, 440.9]
B群（Test）[454.1, 471.6, 464.8, 484.7, 466.4, 458.8, 502.0, 476.6, 481.0, 458.6, 471.8, 481.7, 462.7, 485.6]
Group A: mean=461.34, unbiased sample variance=81.62, n=14 Group B: mean=472.89, unbiased sample variance=174.83, n=14 calculated degree of freedom: 23.0Welch's t_test Results t-value: 2.70 p-value: 0.0128 average diff: 11.55 95%CI: (2.70, 20.41)
統計や機械学習においては、今手元にある観測データは、なんらかの確率分布から生成されたという考え方をする。
例えばこのカフェの場合、ラテ注文者の金額は確率的に決まる（例えば10人中5人くらいは400円の注文をする）という考え方。
例えば実際の注文額（観測データ）は、平均400円、ばらつき50円くらいの正規分布から発生したのだろうといった捉え方ができる。
確率分布とパラメータ
パラメータの選定
古典的統計学の場合
不確実性
パラメータの不確かさを含んだ推定
ベイズ推定の場合
なお尤度はあくまで分析者が仮定した確率モデルのパラメータとの適合度合いを表すだけ。観測データとパラメータとの辻褄がもっとも合う尤度だったからといって、最適解のパラメータであることを常時保証してくれるわけではない。
ベイズ推定のモデリングで分析者が設定すること
観測データの生成プロセス
各パラメータの事前分布
関心のある指標自体もモデリング
# Model Definition
# x_A:Control, x_B:Test
x_A = [465.0, 458.6, 466.5, 475.2, 457.7, 457.7, 475.8, 467.7, 455.3, 465.4, 455.4, 455.3, 462.4, 440.9]
x_B = [454.1, 471.6, 464.8, 484.7, 466.4, 458.8, 502.0, 476.6, 481.0, 458.6, 471.8, 481.7, 462.7, 485.6]
# prior distributions of mu 観測データの全体から設定
combined_data = np.concatenate([x_A, x_B])
pooled_mean = combined_data.mean()
pooled_std = combined_data.std() * 2
# prior distributions of sigma データのスケールに応じて動的に設定
std_low = pooled_std * 0.001
std_high = pooled_std * 10
# prior distributions of nu クルシュケさんのベストプラクティスに準じる
lambda_kruschke = 29
# 平均と標準偏差それぞれの2群間差分、効果量、平均の相対比率についてもモデル作成
with pm.Model() as model_weakly:
mu_A = pm.Normal('mu_A', mu=pooled_mean, sigma=pooled_std)
mu_B = pm.Normal('mu_B', mu=pooled_mean, sigma=pooled_std)
# muが負の値を取るのを完全に防ぐなら、TruncatedNormalを使いlowerを指定する方法もあり
# lower=0よりは、ドメイン知識を反映して、ラテの最低価格にするのがより現実的
# mu_A = pm.TruncatedNormal('mu_A', mu=pooled_mean, sigma=pooled_std, lower=0)
# mu_B = pm.TruncatedNormal('mu_B', mu=pooled_mean, sigma=pooled_std, lower=0)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
# MCMC Sampling
with model_weakly:
prior_weakly = pm.sample_prior_predictive()
trace_weakly = pm.sample(tune=2000, draws=2000, random_seed=42, return_inferencedata=True, chains=4,
# target_accept=0.95, # Uncomment if divergences occur (Raise to 0.9-0.95)
idata_kwargs={"log_likelihood": True}) # モデル比較に使用する対数尤度を保存
posterior_predictive_weakly = pm.sample_posterior_predictive(trace_weakly)
ベイズ推定においては、サンプリング（マルコフ連鎖）が正常に機能し、事後分布に正しく収束しているかを確認する必要がある
trace_plotとsummaryを確認し、基準値をクリアしていれば適切な推定結果と解釈しやすい
trace_plotによる視覚的チェック
summaryによる数値的チェック
plot_autocorrによる自己相関チェック（任意）
# trace_plot
az.plot_trace(trace_weakly)
plt.tight_layout();
# summary
az.summary(trace_weakly)
# mcse_meanが各パラメーターの標準偏差の5%以下（できれば2%以下）という基準値に収まっているかをチェック
summary_df = az.summary(trace_weakly)
summary_df['mcse_mean / sd'] = summary_df['mcse_mean'] / summary_df['sd']
summary_df['mcse_mean / sd (%)'] = (summary_df['mcse_mean'] / summary_df['sd']) * 100
output_df = summary_df[['sd', 'mcse_mean', 'mcse_mean / sd', 'mcse_mean / sd (%)']]
rounded_df = output_df.round({'sd': 1, 'mcse_mean': 3, 'mcse_mean / sd': 4, 'mcse_mean / sd (%)': 2 })
rounded_df
# # 自己相関
# az.plot_autocorr(trace_weakly, var_names=["nu"], figsize=(8,3), textsize=10)
# plt.tight_layout();
# forestplot
# HDIを一覧で視覚的に解釈しやすい（太線は50%HDI）
az.plot_forest(trace_weakly, figsize=(6,4), textsize=10);
読み解き方
本件の場合
az.plot_ppc(posterior_predictive_weakly, num_pp_samples=100, figsize=(10, 4));
# posterior
az.plot_posterior(trace_weakly, figsize=(10,6), textsize=10)
plt.tight_layout();
# posterior
# 差分、効果量、平均の変化率のmeanやhdiの表示
az.plot_posterior(trace_weakly, # hdi_prob=0.95,
var_names=["mean_diff", "relative_mean_diff"],
ref_val=0, color='#87ceeb', figsize=(8,2), textsize=10);
az.plot_posterior(trace_weakly, # hdi_prob=0.95,
var_names=["std_diff", "effect_size"],
ref_val=0, color='#87ceeb', figsize=(8,2), textsize=10);
# 下準備
# 多次元配列を一次元のサンプル配列として抽出
mu_diff_samples = az.extract(trace_weakly, var_names="mean_diff").values
# posterior_meanと94%HDIをグラフに表示するために値取得
posterior_mean = mu_diff_samples.mean()
hdi_lower = np.percentile(mu_diff_samples, 3)
hdi_upper = np.percentile(mu_diff_samples, 97)
# 2群間の平均値の差分の事後分布を表示
fig, ax = plt.subplots(figsize=(6,4))
ax.hist(mu_diff_samples, bins='auto', density=True, # 縦軸を確率密度（全体の面積が1）にする
color="tab:orange", alpha=0.5)
# bins='auto'にすると、データの四分位範囲とサンプルサイズから、外れ値に影響されにくい最適なビン幅を自動計算してくれる
# density=Trueにすると、ヒストグラムの面積の合計が1になる（確率密度関数と同じ状態）になる！
ax.axvline(0, lw=0.8, linestyle="dashed") # ゼロのタテ線
# observed_meanと95%CIの描画（観測データの平均値と95%信頼区間）
# y位置を決める（グラフの一番下から少し上げた高さ）
ymin, ymax = ax.get_ylim()
ci_y = ymin + 0.07 * (ymax - ymin)
ax.scatter(diff_mean, ci_y, s=30, facecolors="white", edgecolors='gray', linewidths=2, label= f"Observed Average: {diff_mean:.2f}" , zorder=5) # mean の点
ax.hlines(y=ci_y, xmin=ci_lower_w, xmax=ci_upper_w, color='gray', linewidth=3, label=f"95% CI: [{ci_lower_w:.2f}, {ci_upper_w:.2f}]") # 95%CIの横線
# posterior_meanと94%HDIの描画（ベイズ推定での事後分布）
# y位置を決める（グラフの一番下から少し上げた高さ）
ymin, ymax = ax.get_ylim()
hdi_y = ymin + 0.03 * (ymax - ymin)
ax.scatter(posterior_mean, hdi_y, s=30, facecolors="white", edgecolors='#7f2704', linewidths=2, label= f"Posterior Average: {posterior_mean:.2f}" , zorder=5) # mean の点
ax.hlines(y=hdi_y, xmin=hdi_lower, xmax=hdi_upper, color='#cc5500', linewidth=3, label=f"94% HDI: [{hdi_lower:.2f}, {hdi_upper:.2f}]") # 94%HDIの横線
ax.legend(loc='upper right')
ax.set_xlabel("Difference in Average (JPY, Test - Control)", fontsize=10)
ax.set_ylabel("Probability Density (Area integrates to 1)", fontsize=10)
ax.set_yticklabels([]) # y軸の目盛りラベルを非表示にする（確率密度関数の値自体に意味を見出さなくていいように）
ax.tick_params(axis='y', which='both', length=0) # y軸の目盛り線の長さを0にして非表示状態にする
# ax.get_yaxis().set_visible(False) # これだとlabelも非表示になる
ax.set_title("Estimated Lift of Latte \n(t_dist Weakly Prior Model)", fontsize=12)
plt.tight_layout()
plt.show()
# print('observed average: {0:,.2f} 95%CI: ({1:,.2f}, {2:,.2f})'.format(diff_mean, ci_lower_w, ci_upper_w))
# print('posterior average: {0:,.2f} 94%HDI: ({1:,.2f}, {2:,.2f})'.format(posterior_mean, hdi_lower, hdi_upper))
# Model Definition
with pm.Model() as model_informative:
# A店は、手元のデータの全体平均をベースに緩やかに置いておく
mu_A = pm.Normal('mu_A', mu=pooled_mean, sigma=pooled_std)
# B店（Test群）の事前分布に「50円プラス」の知識を組み込む。
# 過去の実績（50円アップ）を信じつつ、ブレ（sigma）も考慮する
prior_B_mean = pooled_mean + 50
mu_B = pm.Normal('mu_B', mu=prior_B_mean, sigma=pooled_std)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
with pm.Model() as model_truncated:
# muが負の値を取るのを完全に防ぐなら、TruncatedNormalを使いlowerを指定する方法もあり
# lower=0よりは、ドメイン知識を反映して、ラテの最低価格にするのがより現実的
mu_A = pm.TruncatedNormal('mu_A', mu=pooled_mean, sigma=pooled_std, lower=0)
mu_B = pm.TruncatedNormal('mu_B', mu=pooled_mean, sigma=pooled_std, lower=0)
std_A = pm.Uniform('std_A', lower=std_low, upper=std_high)
std_B = pm.Uniform('std_B', lower=std_low, upper=std_high)
nu_kruschke = pm.Exponential('nu', lam=1/lambda_kruschke) + 1
obs_A = pm.StudentT('obs_A', mu=mu_A, sigma=std_A, nu=nu_kruschke, observed=x_A)
obs_B = pm.StudentT('obs_B', mu=mu_B, sigma=std_B, nu=nu_kruschke, observed=x_B)
mu_diff = pm.Deterministic('mean_diff', mu_B - mu_A)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A)
std_diff = pm.Deterministic('std_diff', std_B - std_A)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A**2 + std_B**2) / 2))
# Model Definition
with pm.Model() as model_gamma:
# ガンマ分布のパラメータは「正の実数」でないとなので弱情報事前分布としてExponentialやHalfNormalを使うことが多い
alpha_A = pm.Exponential('alpha_A', lam=0.1)
beta_A = pm.Exponential('beta_A', lam=0.1)
alpha_B = pm.Exponential('alpha_B', lam=0.1)
beta_B = pm.Exponential('beta_B', lam=0.1)
obs_A = pm.Gamma('obs_A', alpha=alpha_A, beta=beta_A, observed=x_A)
obs_B = pm.Gamma('obs_B', alpha=alpha_B, beta=beta_B, observed=x_B)
# ガンマ分布の数理特性から、事後分布の平均と標準偏差を逆算する（deterministic用）
# 各群の期待値（平均値 = alpha / beta）の計算
mu_A_calc = pm.Deterministic('mu_A_calc', alpha_A / beta_A)
mu_B_calc = pm.Deterministic('mu_B_calc', alpha_B / beta_B)
# 各群の標準偏差（sd = sqrt(alpha) / beta）の計算
std_A_calc = pm.Deterministic('std_A_calc', np.sqrt(alpha_A) / beta_A)
std_B_calc = pm.Deterministic('std_B_calc', np.sqrt(alpha_B) / beta_B)
# 平均と標準偏差それぞれの2群間差分、効果量、平均の相対比率についてもモデル作成
mu_diff = pm.Deterministic('mean_diff', mu_B_calc - mu_A_calc)
mu_diff_rr = pm.Deterministic('relative_mean_diff', mu_diff / mu_A_calc)
std_diff = pm.Deterministic('std_diff', std_B_calc - std_A_calc)
effect_size = pm.Deterministic('effect_size', mu_diff / np.sqrt((std_A_calc**2 + std_B_calc**2) / 2))
# 各モデルのサンプリング結果を格納して、どの観測値の尤度で比較するかを変数名で指定する
model_compare = az.compare({
"t_dist Weakly Prior Model": trace_weakly,
"t_dist Informative Prior Model": trace_informative,
"t_dist Truncated Model": trace_truncated,
"Gamma Weakly Prior Model": trace_gamma
},
var_name="obs_B",
# ic="waic" # デフォルトはloo
)
az.plot_compare(model_compare);
model_compare
『How Brands Grow : What Marketers Don't Know』バイロン・シャープ
マーケティング界では常識のように語られている法則や理論の中には、実はなんの裏付けもなかったり、誇大表現がまかり通っているものがある。
著者バイロン・シャープ Byron Sharpは、従来のマーケティング理論をやみくもに盲信するのではなく、エビデンスに基づいたマーケティングを提唱。実際の販売データや市場調査データの裏付けとともに新たな法則・理論を提示しながら、より効率的なマーケティング活動をおこなうことの重要性を説く。
エビデンスに欠ける従来思想に固執して出血大サービスのマーケティングを続けるのか、それともデータで検証された確からしい見識を携えてブランドを有意義に成長させるのか。見なかったフリをして進化を止めるのか、それとも積極的に取り込んで自らをも成長させるのか。後者でいたい方には有益な書籍。
Evidence-based Marketing
How Brands Grow
How to Grow Your Customer Base
Which Customers Matter Most?
Our Buyers Are Different
Who Do You Really Compete With?
Passionate Consumer Commitment
Differentiation versus Distinctiveness
How Advertising Really Works
What Price Promotions Really Do
Why Loyalty Programs Don't Work
Mental and Physical Availability
A Final Word
ブランドが成長することは可能だ。新しいブランドである必要はない。優れた広告、優れたブランディング、優れたメディア戦略、優れた店内ディスプレイ、そして前述の7つのルールに従うこと、これらすべてが成長への道だ。
和田信明 中田豊一
国際援助においては、現地の人々の合意を得ながら開発プログラムが進む。だが後で振り返ってみると、プログラム内容が最適な解決策になっていなかったり、一過性に終わっている場合もある。その原因の一つに、援助する側とされる側の認識のズレが挙げられる。
例えば「この村では何が問題で困ってますか？」「水汲みが大変です」「井戸はないんですか？」「自分たちでは作れません」「では作るのを手伝いましょう」井戸が作られたおかげで水汲みに遠出しなくて済むようになりました、めでたしめでたし。 ところが数年後に村を訪れると井戸が使われていない。井戸のレバーが壊れたらしい。住民は以前のように水汲みをしている。井戸の修理方法がわからないし、水汲みはもともとやっていたし、そこまで不便とも思っていないらしい。井戸作成プログラムは一定の役には立った。だがさして困っていない様子を見ると、井戸はあれば便利かなぐらいの要望であって、真に解決すべき重大な問題が他にあったのかもしれない。
認識のズレをなくすにはどうすればよいのか？そのカギは相手への質問の仕方にあった。NGOとして途上国に赴いた著者達が、現地の住民やスタッフとの会話で実践しながら編み出したコミュニケーション方法を体系化。実際の出来事を紹介しながら具体的な手順、注意点について解説している。
想起しやすい事実を口にしながら過去の出来事を整理してもらうことで、相手の記憶の棚卸を手伝う。 「その問題が最初に起こったのはいつ？」「最近その問題が発生したのはどこで？」 「その問題で誰が困っている？」「その問題に対してこれまでに対処したことは何？」 事実を重ねる対話を続けることで、相手は気づきを得やすくなる。
「なぜ」「どうして」は相手の負担が大きい、詰問されてるようにも感じる。回答が難しい。 空気を読んで、質問者が欲しがっている答えを口にしてしまう。質問する側も自分の思い込みで質問してしまう。 そうすると事実から遠ざかる。対話の焦点が揃わない。問題が的確に定義されない。解決方法も的外れになる。 事実をたずねることで、本当に解決したい問題の定義、そして的確な答えへとつながる。
書籍タイトルからは遠い途上国の話に見える。だがその応用範囲は広い。日常のあらゆるコミュニケーションに通じる。ビジネスにおけるファシリテーション、問題解決にも使える。
答えを急がない。問題を取り違えない。相手と違う景色を見ない。なぜをたずねるよりも「事実」を着々と積み重ねる。 地味に見えるこの対話術は、途上国に限らずあらゆる場面で有効に機能する。
同著者の手により、実社会において適用しやすいようにエッセンスを抽出した後続書籍あり。おそらくこちらの方が入手しやすく、かつ読みやすいと思われる。
「良い質問」を40年磨き続けた対話のプロがたどり着いた 「なぜ」と聞かない質問術 中田 豊一
照井伸彦
タイトル通り、ベイズ統計モデリングの手法をマーケティング場面に活用することを目指した本
難解な点は、
違う表現をすると、
消費者行動場面を、ベイズ的観点からモデリングしてみる際の手順をちょっと知ってみたい人にはおすすめ 書かれてることを自分の中で咀嚼して応用が効かせられる人にもおすすめ
木戸 茂
マーケティングをデータに基づいて科学的に意思決定する際の手法や応用手順を体系的に紹介するシリーズ書籍。
テレビや紙の媒体でかつて実際におこなわれてきた、以下のような統計的、計量的手法によるモデル構築と検証結果を説明
広告の短期的効果と長期的成果を区別したり、企業ブランドと製品ブランドを区別するなど、シンプルなモデルを目指してるのは好感持てる。広告の継続出稿有意義性を後押しする材料にも使いやすそう
GRPのような広告につきものの説明変数だけでなく、同一市場内での競合との広告シェア、自社の値引率、自社の相対的なシェア、平均購入個数、継続購入率、市場自体の成長率、市場が寡占的か激戦区かなども変数として採用しており、広告出稿に閉じた話に終わらず、ソリューション全体の提案として使いやすい
プリファレンスという1つの説明変数で購買予測するUSJ森岡モデルに比べて、本書はプリファレンスに集約してしまう前の段階、すなわちブランド・エクイティの説明変数が具体的でわかりやすい
またマーケティング全般というよりは広告業界視点で語られている点は、広告業界関係者にはとっつきやすく感じられるだろう
数学的方面からのアプローチや、サイエンス方面からの理解を深めたい人にはおすすめ
統計スポットライト・シリーズ
野外の生き物を観察して得たデータを分析しながら統計手法を学ぶ本。 -クマさんにGPSをつけて計測した移動軌跡をクラスター解析でデータ集約して生態環境を研究
といった事例を紹介しながら、データの計測の仕方、記録方法、データ集約手順や解析を説明。
自然観察から得られる複雑な事象を数値化、図示化するのはなかなか手間に思えますが、そういうのができるソフトの紹介やデータ加工手順については参考になりそう。
具体的なコードなどの記載まではなく、さらっと概要を紹介した薄い本なので、そういうのが欲しい場合は緑本の方が向いてそう。
統計スポットライト・シリーズ
ポアソン分布の理解を深めることを目的にした本。
へーエクセルでそういうこともできるんだ
上の野外調査ともども、緑本と内容がカブってる。緑本は実践に重きを置いてるので原理や理由の詳細説明はあまりない。それを補えると捉えればこの2冊は有用なんだが、副読本の立ち位置でよさそう。
シミュレーションで社会問題を解決しよう、その取り組みや事例を提示、入門となる本。
共有地の悲劇がもたらす環境問題、エージェント・ベース・モデリング、データ分析でシミュレーション、デマツイートとウィルス感染、都市形成、居住民族の構成変化、三次元物体の透視可視化などを紹介しながら、シミュレーションの楽しさや難しさを語る。
いろんなテーマをちょっとずつ紹介してあり、また実生活とも関わり合いの多い事象を取り上げてるので、読みやすかった。
興味を持ってもらうための導入とわりきってもともと作ってあるため、詳しくは参考文献見て自助努力でと放り出される。
放送大学のソーシャルシティ（'17）講座は続きとしてよいかも。 http://www.ouj.ac.jp/hp/kamoku/H30/kyouyou/C/seikatu/1519077.html
放送大学のははシミュレーション関係ない。都市問題つながり。福岡の繁華街、天神に新しい商業施設がオープン、駐車場は足りるの？大分市の駅高架化で、駅ビルがきれいになったら商店街との人の行き来はどう変わった？梅田とか落合氏とかきゅんくんなども登場。
臨床試験での実験計画において従来使われてきた頻度流統計の利点と欠点を述べつつ、それを補う形でベイズ流を説明する本。
頻度流との対比でベイズ流を説明してるので違いやそれぞれのよさがわかりやすい
残念ポイントは特に見当たらない。じゃあNPSが高いかというとそういう気もしない。なんでだろう、これ一冊では完結しないからだろうか。うまく転用することで有効活用できそうな気はする（ちゃんと咀嚼できてない）
Pythonでコレポンをやる
略してコレポン
point
Rの場合
library(ca) ca(table.T)plot(ca(table.T))Pythonの場合
python -m pip install mca理論的には、行・列ともにカテゴリデータの頻度表（人数、出現数などのcountデータ）を対象としている
コレポンをよく使う場面
他に、
比率表（回答者数全体に占める割合、%データ）のデータセット
# import
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import mca
# version
import sys;print("Python", sys.version.split()[0])
import platform;print("OS", platform.platform())
print("is_google_colab:", "google.colab" in sys.modules)
from importlib.metadata import version
packages = [
"pandas",
"numpy",
"matplotlib",
"mca"
]
for pkg in packages:
try:
print(pkg, version(pkg))
except Exception as e:
print(pkg, "not installed")
Python 3.11.15 OS macOS-26.5-arm64-arm-64bit is_google_colab: False pandas 2.2.0 numpy 1.26.4 matplotlib 3.10.9 mca 1.0.4
# サンプルのdataframeを生成する
# メニュー（行）
menus = ['Coffee', 'Tea', 'Latte', 'Espresso', 'Cake', 'Sandwich', 'Cookie', 'Toast']
# イメージワード（列）
words = ['香り高い', 'コスパが良い', '見た目が映える', '濃厚な味わい',
'ボリューム満点', '手軽に食べれる', '贅沢な気分', '健康的なイメージ']
# プロットがある程度バラつくように、『Tea』『Cookie』『手軽に食べれる』が中心近辺にくるように平均的な値に寄せている
data = [
[50, 25, 5, 15, 5, 30, 10, 5], # Coffee
[25, 20, 15, 10, 10, 35, 15, 15], # Tea
[10, 10, 45, 35, 10, 20, 30, 5], # Latte
[40, 10, 10, 55, 5, 15, 25, 5], # Espresso
[10, 10, 50, 30, 10, 15, 45, 5], # Cake
[ 5, 25, 10, 5, 50, 40, 5, 30], # Sandwich
[20, 20, 18, 15, 15, 30, 15, 12], # Cookie
[10, 40, 5, 5, 35, 40, 5, 20], # Toast
]
df = pd.DataFrame(data, index=menus, columns=words)
# MCAの実行
ncol = df.shape[1]
mca_ben = mca.MCA(df, ncols=ncol, benzecri=False, TOL=1e-8)
# 行スコア（=座標）
result_by_row = pd.DataFrame(mca_ben.fs_r(N=2))
result_by_row.index = list(df.index)
print ("Score by row:")
print(result_by_row.round(3))
Score by row:
0 1
Coffee 0.010 -0.686
Tea 0.132 -0.137
Latte -0.513 0.348
Espresso -0.535 -0.413
Cake -0.594 0.440
Sandwich 0.771 0.284
Cookie 0.080 -0.020
Toast 0.710 0.047
# 列スコア（=座標）
result_by_column = pd.DataFrame(mca_ben.fs_c(N=2))
result_by_column.index = list(df.columns)
print ("Score by column:")
print(result_by_column.round(3))
Score by column:
0 1
香り高い -0.184 -0.727
コスパが良い 0.433 -0.130
見た目が映える -0.533 0.537
濃厚な味わい -0.629 -0.123
ボリューム満点 0.724 0.326
手軽に食べれる 0.339 -0.065
贅沢な気分 -0.579 0.227
健康的なイメージ 0.642 0.160
# 固有値と寄与率
# MCAでは「行と列のどちらか少ない方の数-1」 が固有値（eigenvalue）の最大数になるためそれを算出
num_of_eigenvalue = min(len(df.index), len(df.columns)) - 1
data = {
'Factor': range(1, len(mca_ben.L) + 1),
'value': pd.Series(mca_ben.L),
'ratio': mca_ben.expl_var(greenacre=False, N=num_of_eigenvalue)
}
columns = ['Factor', 'value', 'ratio']
table2 = pd.DataFrame(data=data, columns=columns).fillna(0)
# ratioの累計（cumulative_ratio）列を追加
table2['cum_ratio'] = table2['ratio'].cumsum()
print("Principal & Inertia(固有値 寄与率):")
print(table2.round(3).to_string(index=False))
Principal & Inertia(固有値 寄与率):
Factor value ratio cum_ratio
1 0.266 0.598 0.598
2 0.131 0.296 0.894
3 0.033 0.075 0.969
4 0.007 0.016 0.986
5 0.004 0.008 0.994
6 0.003 0.006 1.000
7 0.000 0.000 1.000
# 上記Factor1,2でプロット
plt.figure(figsize=(7, 6))
# plt.rcParams['font.family'] = ['sans-serif', 'Meiryo'] # windows
plt.rcParams['font.family'] = ['sans-serif', 'Hiragino Sans'] # mac
# plot by row
plt.scatter(result_by_row[0], result_by_row[1], color='tab:blue', s=20, marker="o", label="Row")
for k, v in result_by_row.iterrows():
plt.annotate(k, v, xytext=(2, 2), textcoords="offset points")
# plot by column
plt.scatter(result_by_column[0], result_by_column[1], color='tab:orange', s=20, marker='s', facecolors="none", label="Column")
for k, v in result_by_column.iterrows():
plt.annotate(k, v, xytext=(2, 2), textcoords="offset points")
plt.axhline(0, color='gray', linestyle='dashed')
plt.axvline(0, color='gray', linestyle='dashed')
plt.xlabel('Factor 1')
plt.ylabel('Factor 2')
plt.title("Correspondence Analysis\n(Customer Satisfaction Survey Example)")
plt.legend(loc='best')
# plt.gca().set_aspect("equal")
plt.show()
（2026年5月更新）