ラベル プログラミング の投稿を表示しています。すべての投稿を表示

2020年9月30日水曜日

Prophet by Python on Mac


# pip3 install pystan
# pip3 install fbprophet

python3 -m pip install prophet
python3 -m pip install --upgrade plotly

2020年6月14日日曜日

Pythonや機械学習を学ぶ

機械学習の前に重要なデータ抽出・加工に便利なPythonライブラリ「pandas」の基本的な使い方のチュートリアル

https://www.atmarkit.co.jp/ait/articles/1802/13/news012.html

［Python入門］リストの操作 (1/4)

https://www.atmarkit.co.jp/ait/articles/1906/04/news009.html

時系列データへの回帰分析

https://logics-of-blue.com/time-series-regression/

【Python】pandasのDataframe操作

http://canisterism.hatenablog.com/entry/2018/01/07/150105

#特定のセルの値を取り出してseriesに

pythonでexcelファイル処理まとめ

https://qiita.com/hasepy/items/06d5d2e2b6495752442c

[Python] Excelで文字化けしないCSVファイルを書き出す

https://qiita.com/y4m3/items/674423b596284bbc7cf7

Pythonの日本語処理

http://www.wakayama-u.ac.jp/~kazama/lab/python/i18n.html

標準出力の文字コードと自動変換

https://www.javadrive.jp/python/japan/index2.html

文字コードの指定

https://www.javadrive.jp/python/japan/index1.html

Unicode文字列(ユニコード文字列)

https://www.javadrive.jp/python/string/index5.html

日本語文字列コード問題まとめ

http://python.matrix.jp/pages/tips/string/encoding.html

python3で日本語を含むURL(日本語URL)にurllibでアクセスする際、htmlで勝手にencodeされてエラーするので回避策をメモ

https://qiita.com/mix/items/87d094414e46f857de45

日本語を含むURLでつまづく

http://mankuro.xyz/blog/2017/04/25/japanese-url/

単一ドメインを走査する（URLに日本語を含む場合）

http://nnpo.hatenablog.com/entry/2016/12/13/002540

PythonでURLエンコード/デコード

http://shuzo-kino.hateblo.jp/entry/2016/10/22/224034

pythonでutf-8日本語文字列を、URIエンコード・URLデコードする

http://linux.oboe-gaki.com/archives/000333.html

Pythonでマルチバイト文字を扱う際に気をつける点。

https://gist.github.com/devlights/4561968

Python と文字コード

http://www.kabipan.com/computer/python/unicode.html

[Python] urllib.parseによるURLエンコード／デコードの方法

https://hibiki-press.tech/learn_prog/python/url-encoding/2804

機械学習

第1回　難しくない！　PyTorchでニューラルネットワークの基本
https://www.atmarkit.co.jp/ait/articles/2002/06/news025.html

AI・機械学習のための数学超入門 ― 前提知識は四則演算だけ！ (1/4)
https://www.atmarkit.co.jp/ait/articles/2003/02/news023.html

機械学習やディープラーニングってどんなもの？ (1/2)
https://www.atmarkit.co.jp/ait/articles/2003/24/news016.html

「機械学習の最先端」を効率的に情報収集！　おすすめのメルマガ3選
https://www.atmarkit.co.jp/ait/articles/2006/11/news016.html

機械学習の手法13選ー初級者、中級者別に解説！
https://ainow.ai/2020/06/04/222969/

スターバックスはコーヒー事業者ではない ― データテック企業なのだ
https://ainow.ai/2020/05/28/222388/

平均二乗誤差

pytorchにはMSELossがある

criterion = torch.nn.MSELoss()

y = f(x)
loss = criterion(y, t)
print(loss.data)

自作するとこう
def mycriterion(x, y):
result = (x - y) ** 2
return result.sum() / len(result)

myloss = mycriterion(y, t)
print(myloss)

2020年5月24日日曜日

Pythonで体験するベイズ推論

Pythonで体験するベイズ推論 PyMCによるMCMC入門
キャメロンデビッドソン=ピロン（著）　玉木徹（翻訳）

PyMC3のインストール方法

岩波データサイエンス Vol.1
[特集] ベイズ推論とMCMCのフリーソフトのサポートページ
https://sites.google.com/site/iwanamidatascience/vol1/support_tokushu#TOC-PyMC3-

2020年4月12日日曜日

2群の平均値の差の検定そしてベイズ推定

背景

2群の平均値に差はあるか？を簡便におこなうならt検定
中心極限定理とt分布の頑健性のおかげで、たいていの場面において有用
計算簡単でコスパに優れた方法だが、有意差がある/ない以上の情報はもたらしてくれない

そこでベイズ推定

ここでは2群の観測データそれぞれがt分布から生成されると仮定したモデルを作成する
そしてこのモデルにより、2群の平均値の差の事後分布を得る
これにより2群の平均の差が一定値以上である確率など、豊かな情報を得ることができる

結論

ベイズ推定においても、2群の平均の差の平均値や区間推定は、Welchのt検定とほぼ同じ結果が得られた
弱情報事前分布の使用により事後分布はほぼ尤度だけで形作られ、しかも単峰だったため
不確かさを扱いながらの推定は手間を要するが、手間をかけたぶん実務面での問いに答えやすい情報を得られた

シナリオ

背景: カフェチェーンの店頭ディスプレイ改善施策
- あるカフェチェーンで、店内POP（販促物）を刷新することで「ラテ」の売上拡大を企図している
- 立地や客層、売上規模、ラテの販売動向が似通っている2店舗を対象に、ABテストを模した実験をおこなった
データ構造と検証項目
- 2店それぞれの14日間のラテ注文単価データ(日次ラテ売上金額 / 日次ラテ注文数)
- ラテのサイズやトッピングの有無などにより変動しやすい売上金額ではなく注文単価で比較
- 処置（テスト店舗での新POP）は結果変数（ラテの日次注文単価）に影響してそうかを検証したい

観測データ

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]

Welch t-test

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.0

Welch'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円くらいの正規分布から発生したのだろうといった捉え方ができる。
確率分布とパラメータ
- ただし真の確率分布は誰にもわからない（本当は平均380円だったのかも？本当はもっと複雑な形状の分布だったのかも？）。
- 分析者はなるべく真の確率分布に近しそうな確率モデルを構築する。この確率モデルの形状を決めるのがパラメータ。
- 真の確率分布がわからない以上、適切な確率モデルもパラメータも分析者にはわからない。
- 一方、手元には観測データがある。事実としてここにあるからには確からしいであろうと思える。
- そこで、観測データは確からしいであろうを出発点として、このデータを生み出したパラメータを（観測データの生成メカニズムとは）逆向きに探るアプローチをとる。
パラメータの選定
- 手元の観測データを得るのにちょうどいいパラメータもあれば、そうでないものもある。
- 良し悪しの判断に用いるのは、データとパラメータの適合度合いを表す指標（likelihood 尤度）。
- データを固定して、パラメータをいろいろと変化させたときに、尤度がどう変わるかをグラフや数式で表現したのが尤度関数(Likelihood Function)。
古典的統計学の場合
- 尤度が一番しっくりくる箇所を特定してパラメータを推定する（例えば二次関数を微分して極値を求める最尤法など）。t検定の考え方もこれに準ずる。
- 計算が速く推定結果がわかりやすい。しかしそのピンポイントのわかりやすさが仇になる場合もある。
不確実性
- Welchのt検定の結果、新POP導入でTest店のラテ注文単価はControl店を平均的に上回っていたことがわかった。
- だが、14日間のテスト期間中には注文単価が同程度だった日や下回っている日もあった。
- t検定が教えてくれるのは「平均して効果があったかどうか」という1点。
- 14日間という限られた観測データから得られた推定値が、不確実性をどの程度含んでいるかは見えない。
- テスト期間の傾向が今後も続くのか、他店導入時に裏目に出るリスク（単価が下がる確率）はどのくらいかなど実務面での問いには答えにくい。
パラメータの不確かさを含んだ推定
- であれば、今手元にある限られた観測データから生じる「パラメータの不確かさ」を含んだままに推定してみる。
- これにより、まだ見ぬ場面における不確かさを確率的に表現しやすい、使い勝手のいい推定結果を得ることができるのではないか。
- そこで、尤度以外の情報（すでにわかっている知識、過去経験など）も駆使しながら不確かさをうまく扱えるベイズ推定をおこなう。
ベイズ推定の場合
- ベイズ推定では探索によって得られるパラメータを確率分布(事後分布)で表す。
- メリットとして、推定値の不確かさを自然に表現しやすく、実務面での問いにも答えやすい情報を得られる。
1. まずパラメータの探索範囲を確率分布（事前分布）で仮定する。
2. 次に手元の観測データと事前分布を足がかりにして探索する。
3. こうしてパラメータの不確かさを考慮した推定値（事後分布）を得る。
  　※事前分布に、データから得られた尤度（情報）を掛け合わせる（更新する）ことで、事後分布を得る

なお尤度はあくまで分析者が仮定した確率モデルのパラメータとの適合度合いを表すだけ。観測データとパラメータとの辻褄がもっとも合う尤度だったからといって、最適解のパラメータであることを常時保証してくれるわけではない。

Bayesian Model Definition

ベイズ推定のモデリングで分析者が設定すること
- 観測データの生成プロセス（尤度関数となる確率分布の選定）: 観測データがどんなルール（確率モデル）から生まれたか
- 選んだ確率モデルの形状を決めるパラメータの事前分布: 確率モデルの平均やバラツキをコントロールするパラメータの探索範囲
観測データの生成プロセス
- 2群の観測データそれぞれが「t分布」から生成されると想定
  - t分布は、平均mu 標準偏差std 自由度ν（ニュー）の3つのパラメータで形状が決まる
  - t分布は、正規分布よりも両端が厚いぶん、観測データに極端な値があった場合に、たまたま起きたノイズ（標本誤差）として処理しやすい
  - t分布は、νが大きくなるにつれて正規分布に近似していく（t分布を想定することで正規分布をも包含できる）
  - t分布にすることで、小サンプルでも標本誤差を軽減しやすく、不確かさを考慮した「頑健な」推定がしやすい
各パラメータの事前分布
- muとstd共通
  - muとstdの事前分布はかなり広めに設定している（弱情報事前分布 Weakly Informative Prior）
  - 真の値がどのあたりに位置するか確からしい情報がない（観測データを最優先したい、あまり影響与えたくない）との考えを反映
  - この場合、事後分布はほぼ観測データの情報（尤度）だけで形作られる（ので推定結果は検定の結果と近しくなりやすい）
- mu
  - 観測データ全体での平均と、観測データ全体の標準偏差を2倍した標準偏差からなる、裾が長い正規分布としている
- std
  - 観測データ全体のスケールから機械的に算出した、範囲広めの一様分布としている
  - 下限は0付近ギリギリまで探索できるように設定（標準偏差は0以下にはならない、0近辺に張り付く可能性ある）
  - 上限は現実的範囲に制限
- ν（ニュー）
  - 頑健な推定のためにνの値が小さい（＝両端が厚い分布になる）方に高い事前確率を割り当てたい
  - 0をピークに右に裾が長い指数分布は好都合だが、νは0より大きい必要がある。またt分布はνが30を超えるとほぼ正規分布になる（ν=1はCauchy分布と同じ）
  - 以上から平均29（ラムダλ=1/29）の指数分布を右に1だけシフトしている
  - これにより観測データの外れ値を許容しやすいν=1近辺（典型的なt分布）から、ν=30近辺（ほぼ正規分布）まで広く探索できる
関心のある指標自体もモデリング
- 観測データの生成プロセスを推定するついでに、結果の判断や意思決定に有用な情報もモデル化する
- ここでは以下4指標をDeterministic (各種パラメータから導出される決定論的変数)に推定する
  - 平均の差分: リフト差分（Test群の平均 - Control群の平均）
  - 標準偏差の差分: バラツキの差（Test群の標準偏差 - Control群の標準偏差）
  - 効果量: 差の大きさがバラツキに対してどれだけインパクト（大きな効果、意味を持つ差）があるか
  - 平均の変化率: リフトの比率（Control群を基準としたときの相対比率 1=100%すなわち等倍）

# 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)

Sampling Quality Check

ベイズ推定においては、サンプリング（マルコフ連鎖）が正常に機能し、事後分布に正しく収束しているかを確認する必要がある
trace_plotとsummaryを確認し、基準値をクリアしていれば適切な推定結果と解釈しやすい
trace_plotによる視覚的チェック
- 左側はKDE（カーネル密度推定）による事後分布の形状、右側はサンプリングの軌跡（トレースライン）。いずれのグラフもchainsで指定した本数（例: 4本）の線がプロットされる
- 左側のグラフでは4本の線がだいたい同じ形に、右側のグラフでは毛虫のような見た目（縦軸の値が一定範囲を何度も行き来して変動幅がだいたい同じくらい）になっていれば概ね問題ない
- もしも左右グラフの下部に縦のバーが表示されている場合は、divergence 発散が起きていることを表す（縦バーの箇所ではサンプラーが空間の急斜面で足を踏み外してる、シミュレーションが真の軌道から大きく外れてる、つまりサンプリングがうまくいっていないことを表す）
  - 対処策1: target_acceptの値を0.90〜0.95程度に上げてみる（歩幅を細かくして慎重に探索させる）
  - 対処策2: tuneの回数を増やしてみる（ウォーミングアップ期間を伸ばして、地形の学習を確実にする）
summaryによる数値的チェック
- r_hat
  - chain間の不一致度を表す
  - 目安として1.05以下（理想は1.00）であれば問題ない
  - そうなってない=各chainが同じ分布にたどり着いていない（収束していない）、つまりサンプリングがうまくいっていない
- mcse_mean, mcse_sd（マルコフ連鎖モンテカルロ誤差）
  - サンプリングのブレが原因で生じる推定値の誤差（の平均と標準偏差）を表す
  - 目安としてmcse_meanが0.01以下（パラメータの変動範囲が0〜1の場合）、またはパラメータの標準偏差の5%以下（できれば2%以下）と十分に小さければ問題ない
- ess_bulk, ess_tail（有効サンプルサイズ）
  - 自己相関を考慮した実質的なサンプル数
  - 事後分布の形状やHDIを安定して推定できているかを表す
  - 目安として400以上（できれば数千）確保されていれば問題ない
plot_autocorrによる自己相関チェック（任意）
- summaryでess_bulkとr_hatが基準値をクリアしていれば自己相関は心配無用なので実行せずともよい
- 問題ありの場合は、該当のパラメータのみを指定して実行（することでグラフが少ないぶん見やすい）
- どこに問題あるかの原因の切り分けに役立つ
  - 特定のchainだけ自己相関が高い→tuneを増やす、サンプリングの初期化オプションを調整するなど
  - どのchainも自己相関が高い→モデルの構造を変える（パラメータ間の相関などが疑われる）

# 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);

PPC

Posterior Predictive Check: 作成したモデルが観測データをどの程度再現できているかを確認

事後分布のchainから任意のパラメータをランダムに取り出す（例: mu_A=450 std_A=20 nu=5）
取り出したパラメータから成る分布から、観測データと同じサンプルサイズの擬似データを生成する（t分布から14サンプル取り出す）
このKDE（カーネル密度推定）をプロットしたのが青線
ここまでを指定の数（num_pp_samples=100）繰り返す。観測データのサンプルサイズが多い場合や分布次第では、青線が帯のように密集状態になる
青線の平均をとったのがオレンジ点線（=モデルの平均的な予測）

読み解き方
- 黒線（観測データのKDE）が青線の変動範囲やオレンジ線と同じような位置、形であれば、モデルは観測データをよく近似していると解釈できる
- 大きくズレていたら、分布の選定や事前分布の見直し、例えばt分布以外がよかったかな？事前分布の探索範囲を変えた方がよい？など、モデル修正方針の検討材料になる
- あくまで今手元にある観測データの発生メカニズムを矛盾なくシミュレーションできているかを見てるだけ。観測データやモデル自体の確からしさを裏付けるものではない（観測データが偏ってる可能性や過学習リスクは残る）
本件の場合
- 観測データのサンプルサイズが少ないため、青線はあまり密集していない
- 黒線がモデルの予測範囲（青線たちの広がり）の中に収まっており、いい感じのモデルといえそう
- ダメなモデルだったら、黒線が青線たちの外側に大きく飛び出していたり、山の位置が全くズレてたりするだろう
- muの事前分布をNormalにしたことで、ごく僅かに左端がマイナスの値を取る。これを100サンプルが引き当てた場合には、違和感のあるプロットになる。2群の平均値の差の事後分布にはほぼ影響ないが、どうにも気になる場合はTruncatedモデルを採用する手もある

az.plot_ppc(posterior_predictive_weakly, num_pp_samples=100, figsize=(10, 4));

Interpretation

Posterior Inference and Interpretation: 推論された事後分布の解釈
plot_posterior
- 事後分布のKDE, mean,HDIを一覧表示。視覚的に判断しやすい
- ref_valで指定した値を基準に、その値を上回る（下回る）確率を表示することも可能
2群の平均値の差について詳しく
- もっとも関心のあるところなので加工してグラフにいろいろ装飾している

# 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))

Informative Model

ベイズ推定っぽく、事前分布に事前情報を適用したモデル
他チェーンでは新POPで単価が50円上がったらしいという事前情報を反映した事前分布に書き換える
「Test店はControl店よりもあらかじめ50円くらい高くなるポテンシャルを秘めているはずだ」との考えを反映したモデル
尤度関数は引き続きt分布を想定
想定される弱事前情報分布モデルからの変化
- Test店の平均値と94%HDIが上ブレする。ベイズ推定は「他チェーンでは50円上がったらしいけど、今回のデータは11円か…じゃあ間をとってこれくらいが真の実力かな」という、データの不確実性と過去の信頼性を天秤にかけた、いわゆる「収縮効果（Shrinkage）」を起こす。
- ひょっとすると94%HDIが狭まるかも。ベイズ推定は「情報が増えて、より確信を持って（不確かさを減らして）真の値を予測できるぞ」という状態になるので、事後分布の山が尖って、94%HDIの横棒が95%CIに比べて狭く引き締まる可能性がある。

# 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))

Truncated Model

結果変数が負の値をとらないようにするモデル
弱情報事前分布モデルも事前情報を組み込んだモデルも、muの事前分布をNormalにしていた
推定への影響は軽微だが、現実世界のドメイン知識をより厳密に反映したモデルにする
尤度関数は引き続きt分布を想定しつつ、muの事前分布をTruncated Normalに変更する
ドメイン知識を反映したモデルで納得感高めることを企図（副次的にPPCでの違和感ある出力回避可能）

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))

Gamma Model

Truncated Model同様に、結果変数が負の値をとらないようにするモデル
尤度関数はガンマ分布を想定
負の値をとらない連続値であれば、t分布と切断分布の組み合わせよりは、こちらの方が汎用性は高いかも

# 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))

Compare Models

どの尤度関数（モデル）が優れているかを比較
事後分布を丸ごと使って比較する
LOOやWAICといったベイズ専用の指標を計算する。
単に手元のデータへのフィット度を見るだけでなく、「パラメータの不確実性も含めた上で、未来の予測能力が一番高いのはどのモデルか？」を厳密にスコアリングしてくれる。
基本的には WAIC よりも LOO を使う方が推奨とされている。
- 理由は数理的な頑健性（ロバストさ）にある。
- WAIC: 計算がめちゃくちゃ高速だけど、データが少なかったり、極端な外れ値があったりすると、ペナルティ項（p_waic）の計算が少し不安定になって、スコアがブレやすいという弱点がある。
- LOO（正確にはPSIS-LOO）: WAIC とほぼ同じ予測力を測る指標なんだけど、内部で「外れ値に対する警告システム（Pareto k 統計量）」が動いてくれるから、14サンプルみたいな少サンプルデータに対して、より安全で正確に機能する。

# 各モデルのサンプリング結果を格納して、どの観測値の尤度で比較するかを変数名で指定する
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

表の見方

rank: どのモデルが優秀か。
- 予測性能がいい順番。
elpd_loo(elpd_waic): 予測性能のトータルスコア。
- 期待対数予測密度（対数予測密度の期待値 Expected Log Pointwise Predictive Density）。値が大きい、つまり0に近いほど予測性能が優秀。未来のデータをどれだけ精度よく予測できるかを表す。
- t分布たちはほぼ互角だけど、ガンマ分布はかなり低いスコア。つまり「14サンプルの手元のデータに対しては、t分布のほうが圧倒的に座りが良かった」ということを意味している
p_loo(p_waic): 有効パラメータ数。
- モデルの複雑さ（過学習のしやすさ、ペナルティの大きさ）を表す。
- パラメータ数が多かったり、観測データに無理やり形を合わせよう（過剰適合）とすると、この数値が大きくなりelpdが減点される仕組み。
- t分布勢に対してガンマ分布は小さい。ガンマ分布モデルは無駄な複雑さがなく、シンプルにデータを表現していることを示している
elpd_diff: 1位との予測スコアの差分。
- 差が4以上あると有意な差とみなされることが多い
- t分布たちの差分は微量なので1位とほぼ同等と言える。ガンマ分布はずいぶん離れており、今回データにおいては予測力不足と言える。
weight: 予測ブレンドの重み付け
- もしこれらのモデルを合体させて未来を予測するなら、どのモデルの意見を何割信じるべき？というブレンド比率（0〜1）。
- LOOの場合（内部ではスタッキング法という技術が使われることが多い）、モデルの実力が伯仲していたときに0.6と0.4みたいに綺麗に分散して出力されやすい（複数モデルの予測パワーのバランスをリアルに反映しやすい。重み付けの計算がより安定的になる）
- WAICの場合は微差であっても1位のモデルに 1.0（100%）と極端に全振りされやすい。
se: トータルスコアの標準誤差
- 各モデルのelpd_loo(elpd_waic)が、サンプリングのブレによってどれくらいバラつく可能性があるかを表す不確実性の幅（標準誤差）。
dse: 差分の標準誤差（elpd_diff自体の標準誤差）。
- elpd_diff ÷ dseの計算結果が2を超えているかどうかが、実務上の超重要ラインになる。
- 2未満であれば、1位のモデルと数理的には完全に互角・誤差の範囲内。
- ガンマ分布は2を大きく超えており、ガンマ分布が負けているのはサンプリングの偶然ではなく明確に実力差であると言い切れる。
warning: 計算が破綻していないかの警告
- データの中にLOO(WAIC)の前提を揺るがすような過度なばらつき（極端な値）があり、指標の計算が数理的に変になってないかのフラグ。
- True だと数値が信用できない、False だと健康に計算できたという意味。
- LOOのwarning=Trueは、
  - 内部でパレートk統計量（Pareto k diagnostics）という厳密なデバッグシステムが動いて警告を出す。
  - データ1件1件に対して、このデータはモデルにとって予測しづらい外れ値になってないか？を全件チェックして、基準値（k>0.7）を超えたデータが1件でもあるとwarning=Trueになる。
- WAICのwarningは、
  - データ全体のバラつきをなんとなく見て、なんか破綻してそうってアバウトに警告を出すだけ。
  - WAICではt分布勢はサンプルの少なさや裾の厚さのせいでTrueだが、ガンマ分布モデルは Falseだった。
  - ガンマ分布は非負のデータ構造を数理的に正しく捉えているから、警告を出さずにものすごく健康に計算できたという、ガンマ分布の構造的な美しさが証明されている。
scale: 計算の単位空間。
- 対数尤度をどのスケールで計算したか。ArviZでは一律でlog（対数空間）が使われる

モデル評価

4モデルの中ではTruncated Modelが最も予測能力が高いと評価された（僅差のためサンプリング次第でTOP3の並びは変わる可能性あり）
t分布の3モデルはほぼ同程度に対し、ガンマ分布モデルは予測能力が劣ると評価された
ガンマ分布モデルは微量ながらdivergenceが発生していた。また94%HDIの幅が広すぎる推定だった
以上からすると、ガンマ分布モデルを採用するのはなさそう。t分布についてはどれを採用してもよさそうとの結論

2019年1月5日土曜日

画像の圧縮

画像の大きさはそのままに、ファイルサイズを小さくする

ImageMagickでできる

$mogrify -quality 85 *.jpg

ただし一度プレビューで品質下げて書き出すなどしたのを改めてやると、ファイルサイズが大きくなることもある

コマンドで画像を一括処理！ ImageMagickの便利な使い方
https://liginc.co.jp/394506

画像を最適化する

https://developers.google.com/speed/docs/insights/OptimizeImages?hl=ja

フォルダ内のjpegファイルをjpgにする
for file in *.jpeg ; do mv $file ${file/%.jpeg/.jpg} ; done

2018年12月29日土曜日

PythonでMCA（コレスポンデンス分析）

Pythonでコレポンをやる

コレスポンデンス分析

略してコレポン
- 数量化III類（Quantification Theory III）、対応分析（CA: Correspondence Analysis）、多重対応分析（MCA: Multiple Correspondence Analysis）、双対尺度法（Dual Scaling）など、呼び方や近縁手法が複数存在し、数理的には非常に近い関係にある
- ここでは「コレポン」表現をメインに据えている。実際にはCAだったりMCAだったりするが数理的背景は非常に近いため細かい区別は気にしない
point
- 観測データを単純化することで、データの背後にある潜在的な要素を解釈しやすくする
- 質的データ（カテゴリデータ）の関係構造を探索するのに適している
- 主成分分析（Principal Component Analysis）と数理的に近い構造を持つ
- PCA同様にデータを縮約して、2次元プロットなどで可視化できる
- 統計的厳密性よりは、結果解釈のしやすさや仮説探索的手法として喜ばれやすい（コレポンに限らない）
- 大雑把な傾向把握にとどめるのが望ましい（細かいことを気にしないことこそコレポンらしさかも）
Rの場合
1. インストールして library(ca)
2. 読み込んで実行して ca(table.T)
3. 可視化 plot(ca(table.T))
Pythonの場合
- ここではmcaライブラリを使用 python -m pip install mca
- 単機能でシンプルではあるが、新規に始めるならprinceライブラリの方がいろんなことに使いやすいかも

コレポンとPCA

コレポンの特徴
- カテゴリデータや頻度データの関連構造を少数の軸で表現する手法。期待頻度からのズレ（カイ二乗統計量）が大きい方向を探す。
- カイ二乗距離に基づく幾何を前提としている。行・列の周辺度数で基準化した行列を特異値分解する。
- 行と列を双対的に扱い、同一空間上に同時布置できる。
PCA（主成分分析）の特徴
- 数値データの分散構造を少数の軸で表現する手法。データのばらつき（分散）が最大になる方向を探す。
- ユークリッド距離に基づく幾何を前提としている。共分散行列または相関行列を固有値分解する。
- 行（サンプル）は主成分スコア、列（変数、特徴量）は負荷量として表現される。
コレポンとPCAの共通点
- 多変量データの構造を少数の軸で要約・可視化する次元圧縮手法（データの背後にある構造を捉えるためのもの）
- データ中の「似たパターン」や「関係構造」を低次元空間へ埋め込む（低次元に圧縮して、似た回答傾向を可視化する）
- 元データを座標空間へ変換し、近い点ほど似た性質を持つと解釈できる
- 軸ごとの構造を使って、クラスタ傾向や潜在的な特徴を読み解ける
コレポンとPCAの相違点
- コレポンとPCAは少なくとも数理的な骨格は同じなんだが
- 軸の解釈に対する考え方の違い、解析結果の安定性など運用上は多くの相違点が存在する（PCAにおける軸は量的な増減方向として解釈しやすい。コレポンにおける軸はカテゴリの対立軸として表れやすい）
- だがコレポンは適用範囲が広い。理論的前提から多少外れていても、大まかな構造把握用途では有用な場合が多い。特に順序尺度のアンケートデータなどでは、PCA的に見るかコレポン的に見るかの境界が曖昧になることも多い

コレポンの注意点

仮説探索や構造把握に向いた手法であり、推定や検定のような厳密性を主目的とする分析とは少し性格が異なる
適用範囲が広い（深く捉えなくともなんでも実行しやすい）とはいえ、本来持っていた情報を相当捨てていることは気に留めておく
変数間で尺度の意味が大きく異なるデータには向かない（例えば年収と身長と気温を同列の変数として扱うのは不適切）
PCA同様に異常値や希少カテゴリは空間構造をゆがめやすい

コレポンの実情

特にマーケティング分野におけるアンケート調査データにおいて重宝されやすい
「結果の見た目がわかりやすい」「直感的に理解しやすい」「解釈やストーリーを考察しやすいので話が弾む」などが主な理由として考えられる
他に「クロス数表があればできちゃう」「たいていの設問形式に対応できちゃう（ように感じやすい）」「簡単な操作で簡単なアウトプットが得られちゃう」もありそう
理論的美しさよりも、実行しやすさや適用範囲の広さ、データの背後構造の見えやすさが魅力と言える

コレポンの守備範囲

理論的には、行・列ともにカテゴリデータの頻度表（人数、出現数などのcountデータ）を対象としている
- 例えば、行に「カフェのメニュー」、列に「性年代区分」を置いたデータセット（各メニューを好きと回答した人数の性年代別集計値）でおこなうのが王道
- だがクロス数表を読み慣れている場合はコレポンやる意義を感じにくいかも
- 性別や年代による違いは想像がしやすいため（ある程度わかってる顕在的要素をなぞってるだけに見えるぶんコレポンの利点が活かされにくい）
コレポンをよく使う場面
- 例えば、行に「カフェのメニュー」、列に「イメージワード」を置いたデータセット（各メニューの印象にあてはまるワードを選択した人数のワード別集計値）
- イメージが似てるメニューはどれとどれなのか、逆に似ていないのはどれかといった解釈がしやすい（クロス数表からは読み解きにくい）
- コレポンの理論的枠組みからはやや外れるが、コレポンやる意義を感じさせやすい
他に、
- 行に「回答者」、列に「イメージワード」を置いたデータセット（回答者各自がカフェの印象にあてはまるワードに5段階評価したローデータ）
- これも理論的枠組みからは外れる（PCAっぽいともいえる）が、ある人がどんなワードを重視してるか、似た者同士は誰かといった解釈がしやすい
比率表（回答者数全体に占める割合、%データ）のデータセット
- コレポン自体は実行できるし、頻度表よりも入手しやすい場面もあるなどこちらでおこなう人が多数派かもしれない
- これまた理論的枠組みからは外れるが、単純化するという性質からさほど弊害が表れにくい（頻度表での実行結果と大差ないなど）
- 同様理由で、表の縦%や横%が100を超える場合（複数回答形式設問）でもおこなわれているだろう

コレポン結論

理論を超えて使いやすいが、わかりやすさに引きづられない
簡単に得られる情報と引き換えに、捨てている情報も多いことを忘れない

コレポン実行例

あるカフェのお客様アンケート結果をコレポンしてみる
メニュー × イメージワード
- 実務においてありがちなコレポンの使い方
- お店のメニュー（8品）それぞれに対して、印象にあてはまるイメージワード（8項目）をMA（複数回答可能形式）で聴取
- 行（メニュー）、列（イメージワード）のクロス数表（選択した人数の頻度表）
- MA選択率よりはイメージワード各5段階評価のTOP2人数の方がコレポン理論には近い。とはいえ回答者負荷の高い設問形式を推奨したくない意図

# 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()

結果の解釈例

左上のCakeとLatteは似たイメージでお客様から捉えられている
CakeとLatteはともに「映える」「贅沢」といったリッチなイメージが強そう
これらとは対極の右側にはSandwich, Toastが「手軽に食べれる」日常的なMenuとの評価
横軸（第1軸）はハレと日常を表すような背後要因が考えられる
下方のCoffeeはハレと日常の中間的ポジションながらも、Espressoともども「香り高い」イメージ
縦軸（第2軸）は寄与率が1軸ほど大きくなく、背後要因はやや読み解きにくい（今回使用したイメージワードでは不足だった可能性も考えられる）
中心付近にプロットされた「Cookie」「Tea」は他のMenuほどには強い特徴がみられなかったと解釈できる

（2026年5月更新）

2018年11月17日土曜日

exif情報の一括削除

ImageMagickでできる

フォルダ内のExif情報を削除したい場合はオプションの「-strip」で可能

該当フォルダに移動した後、一括操作の「mogrify」を使う

mogrify -strip *.jpg

イメージマジックでフォルダ内のExif情報を一括削除する

https://www.s-oj.com/ec-business/system/imagemagick%E3%82%92%E4%BD%BF%E3%81%A3%E3%81%A6%E3%82%B5%E3%83%96%E3%83%87%E3%82%A3%E3%83%AC%E3%82%AF%E3%83%88%E3%83%AA%E3%81%AEexif%E6%83%85%E5%A0%B1%E3%82%92%E4%B8%80%E6%8B%AC%E5%89%8A%E9%99%A4/

----

すべてのpng画像をjpgに変更

mogrify -format jpg *.png

（元画像は残る）

https://liginc.co.jp/394506

----

ImageMagick基本編：Mogrifyを用いた画像の一括処理について

http://shumilinux.blogspot.com/2016/06/imagemagick-mogrify.html

2017年12月28日木曜日

MacBookAirでPython3環境設定からのexif撮影日時で画像ファイルrename

たぶんここらへん見て環境設定した。

http://qiita.com/yoshizaki_kkgk/items/4663148a2b3ca078ddbc

決定版】MacでPythonを使って『機械学習』を学ぶための環境構築

http://www.sejuku.net/blog/curriculums-python

Python入門完全攻略ガイド【基礎学習/アプリ開発/仕事獲得】

http://qiita.com/yhashizaki/items/c996297a80024bcaa2eb

Pythonスクレイピング（Beautiful Soup）+XPathでテキスト文指定でリンク先ＵＲＬを取得

【参考にさせていただいたサイト】

http://gci.t.u-tokyo.ac.jp/tutorial/crawling/

ウェブページからのデータ取得

http://www.slideshare.net/tushuhei/python-xpath

Python と Xpath でウェブからデータをあつめる

---------------

exiv2 -Pk hoge.jpg
exif情報の種類が見れる

magick identify -verbose hoge.jpg

詳細情報が見れる

---------------

写真アプリの画像ファイルを整理したかった。

撮影日時をexifから取り出してファイル名につけたかった。

ImageMagickというのでできるっぽい

それをinstallしたい

Mac OS X El Capitan - HomebrewでImageMagickのインストール

https://qiita.com/noir/items/dc7366be57f987c04a20

brew install imagemagick

とやったら、command line toolsがないとだめって言われたので

Xcodeをインストールせずに、Command Line Toolsをインストールする方法

https://qiita.com/iwaseasahi/items/eb820762600c815ab100

---------------

iPhotoで、撮影日時をタイトル名に反映させて、タイトル名で書き出すという手もあるが。

https://origin-discussions-jp.apple.com/thread/10183849

---------------

autometerでやる場合

写真のファイル名を撮影した日時にする方法【Mac】

http://www.umurausu.info/blog/archives/add_date_time_to_photo_filename.html

===============

bash 画像の撮影日時（EXIF）を「YYYYMMDD_HHMMSS」形式で取り出してファイル名に付ける（awkとsedで）

http://min117.hatenablog.com/entry/2017/03/04/134130

【ImageMagick】identifyコマンドが実行できない

https://teratail.com/questions/92619

magick identify -verbose IMG_hogehoge.JPG | grep DateTime: | awk '{print $2"_"$3}' | sed -e s/://g

エラーメッセージが出る

identify: unable to open image ' ': No such file or directory @ error/blob.c/OpenBlob/3323.

identify: no decode delegate for this image format `' @ error/constitute.c/ReadImage/509.

20091124_085544

でもexif情報は取り出せてるよね

-------------------------

いろいろやって

それっぽいのを作ってみたが

./rename.sh

が

Permission denied

になるときは

chmod +x rename.sh

とするらしい

これはファイルに実行権限がついていないため、起きる現象です。

chmodコマンドを使って実行権限をつけてやりましょう。

http://wadap.hatenablog.com/entry/20080223/1203787519

rename.shの中身

for i in IMG_*jpg

do if

exiv2 -Pk "$i" | grep -q Exif.Photo.DateTimeOriginal ; then

exiv2 -Pk "$i" | grep -q Exif.Photo.DateTimeOriginal |

~~exiv2 -r'%Y%m%d_%H%M%S'"_:basename:" "$i"~~

exiv2 -r'%Y%m%d'"_:basename:" "$i"

done

んで

./rename.sh

できた！

Exiv2 exception in print action for file IMG_****.jpg:

Failed to read image data

みたいなメッセージが終了時に表示されるがまあいいや

---------------------

magick identify -verbose IMG_****.jpg | grep DateTime: | awk '{print $2"_"$3}' | sed -e s/://g

awk 列を抜き出す
例: 2列目だけ
$ cat << __EOF__ | awk '{print $2}'
1 hoge
2 hage
__EOF__
以下出力
hoge
hage

■ ファイル内の文字列をまとめて置換する
カレントディレクトリのファイル全てを対象にする場合
$ grep -l '置換対象の文字列' 置換対象のファイル | xargs sed -i.bak -e 's/置換対象の文字列/置換後の文字列/g'
https://qiita.com/kkyouhei/items/b4ff839a2f36ba194df3

--------------------

http://microgroove.jp/shaolin/2010/07/osx_automator_script_for_renaming_files_using_exif_datetime.html

http://min117.hatenablog.com/entry/2017/03/04/134130

http://serial-experiments.jp/blog/?p=512

https://blanktar.jp/blog/2015/09/linux-photo-change-datetime.html

http://uxmilk.jp/8366

https://qiita.com/kkyouhei/items/b4ff839a2f36ba194df3

http://www.sw-mono.blog/entry/2017/07/30/115350

http://www.nemotos.net/?p=864

http://www.nemotos.net/?p=745

http://www.nemotos.net/?p=182

http://shellscript.sunone.me/for.html

写真のファイル名を撮影した日時にする方法【Mac】
http://www.umurausu.info/blog/archives/add_date_time_to_photo_filename.html

写真取り込み＆EXIF情報にあわせてファイル名を変更
http://microgroove.jp/shaolin/2010/07/osx_automator_script_for_renaming_files_using_exif_datetime.html

【覚書】Macで写真のリネームをする
http://zipmouse.hatenablog.com/entry/2017/10/07/202640

Mac OS X Mountain Lion でファイルのタイムスタンプ(日付)を変更する方法
http://blog.dreamhive.co.jp/yama/9708.html

bash 画像の撮影日時（EXIF）を「YYYYMMDD_HHMMSS」形式で取り出してファイル名に付ける（awkとsedで）
http://min117.hatenablog.com/entry/2017/03/04/134130

登録: 投稿 (Atom)

ページ

2020年9月30日水曜日

2020年6月14日日曜日

機械学習

平均二乗誤差

2020年5月24日日曜日

PyMC3のインストール方法

Python ヒッチハイク・ガイド

PcMC3のコード

Chapter1の一部をやってみた

2020年4月12日日曜日

シナリオ

観測データ

Welch t-test

統計モデリングのスタンス

Bayesian Model Definition

Sampling Quality Check

PPC

Interpretation

Informative Model

Truncated Model

Gamma Model

Compare Models

表の見方

モデル評価

2019年1月5日土曜日

2018年12月29日土曜日

コレスポンデンス分析

コレポンとPCA

コレポンの注意点

コレポンの実情

コレポンの守備範囲

コレポン結論

コレポン実行例

結果の解釈例

2018年11月17日土曜日

2017年12月28日木曜日