2025-09-01

「決意を新たにする」のはやめよう

ポエム

大事なのは「続けること」よりも「辞めないこと」、の続き。

blog.mudatobunka.org

「決意を新たにする」ってやつは罠なのでしないようにしよう。

「決意を新たにする」という脳内完結型イベント

何が罠かというと、決意したところで "現実の物理的世界は何も変わらない" ことだ。
一方で内面的には、何らかの達成感と満足感が得られてしまう。それで満足してしまえば、ただ脳内で決意して満足して終わっただけになってしまう。

多くの人は達成感が得られていないからこそ・満足していないからこそ手を動かすわけで。まだ何もしてない段階で達成感や満足感が得られてしまう「決意」ってやつはヤバい。

もちろん決意を新たにして行動を起こし、持続すれば素晴らしい。でも行動を起こすのに決意って必須でもないのでは。部屋を片付けたかったら「まず5分だけやるつもりで着手してみろ」ってやつ。

決意だけして満足してしまうリスクを冒すより、もっと軽率にもっと漫然と始めた方がいい。

「辞めないこと」は軽率に再開すること

前回の投稿でも書いたが。「辞めないこと」のポイントは「いつだって軽率に再開していい」というところにある。

物事が続かなかったとき、別に「挫折」というレッテルを貼る必要はない。あくまでも中断してただけだし軽率に再開しただけ、そういう認識で大丈夫。
そこには挫折ないし、再開するためには新たな決意も必要ない。

まとめ

職場で全員に決意表明させるとかも害があるのでは？と思う。自分の決意を言語化してみんなの前で表明する、それだけでエネルギー消費しちゃうし。それで満足しちゃったら、いざ走り出すときに力出ないって。

私からは以上です。

2025-08-29

ベイジアン鶴亀算/t分布モデル

Stan 統計

鶴亀算を解こう。

ある庭に何羽かの鶴と何匹かの亀がいる。
20人の観測者が頭の数と足の数を数えたところ次のようなデータが得られた。

observed_heads observed_feet

96 275

96 272

94 274

95 272

...中略... ...中略...

96 273

庭にいる鶴の数と亀の数を求めよ。

観測者によって頭の数と足の数がバラバラだ。これは各観測者が数え間違いをしている (かもしれない) からだ。

ベイズの考えを使って解く

典型的な鶴亀算の問題は、与えられた「頭の数」と「足の数」を真の値として決定論的に計算する。
今回の問題では「頭の数」と「足の数」に観測誤差を含むので、誤差を加味して「鶴の数」と「亀の数」を評価しなければいけない。

このような問題にはベイズ推定の考え方が使える。隠れたパラメータがあり、パラメータに応じて確率的にデータが観測されていると考える。ベイズの定理 (逆確率の定理) を使うと、データから逆算してパラメータを推定できる。

いま「鶴の数」「亀の数」は未知だが、庭にいる「鶴の数」「亀の数」に応じて「頭の数」「足の数」が観測されているはずだ。ベイズ推定で、観測された「頭の数」「足の数」から「鶴の数」「亀の数」を推定できる。

実行環境・ツール

今回も Stan を使って推定する。Stan は Python から呼び出すことにする。Python の実行環境には Google Colab を使う。
というわけで完全に動作する Google Colab Notebook がここにある。

colab.research.google.com

モデル定義

確率モデルを定義しよう。
ここでの確率モデルとは、「頭の数」「足の数」がどのような分布から生成されているかを仮定するものだと考えればいい。

真の鶴の数 T* と真の亀の数 K* に対して、真の頭の数 H*, 真の足の数 F* は、

H* = T* + K*
F* = 2 T* + 4 K*

こうなるはずだ。

よって観測された頭の数 H と観測された足の数 F は、それぞれ T* + K* と 2 T* + 4 K* を中心とした山なりの分布から生成されていると考えたい。
山なりの分布なので中心である T* + K* や 2 T* + 4 K* の値が観測されやすい。しかし確率分布からの標本なので、もちろん中心から外れてしまうこともある。

今回は H, F が t分布に従うと仮定する。
t分布は3つのパラメータ (母数) によって形が決まる。 $\nu$ , $\mu$ , $\sigma$ だ。

$H \sim StudentT(h | \nu = nu, \mu = T + K, \sigma = sigma\_H)$
$F \sim StudentT(h | \nu = nu, \mu = 2T + 4K, \sigma = sigma\_F)$

このような仮定をおいて T, K, sigma_H, sigma_F, nu を推定対象としよう。

Stan コード

モデルを Stan のコードで書き下すとこうなる。
このコードでは H = T + K, F = 2 T + 4 K という知識を反映しているものの。T = ~, K = ~ の形に解いた式は入っていない。
H, F の値が尤もらしくなるように推定すれば、自然に T, K の尤もらしい値も求まるはずだ。

data {
  int<lower=1> N;      // データ数
  vector[N] H;         // 観測された頭の数
  vector[N] F;         // 観測された足の数
  real H_bar;          // 頭の数の標本平均
}
parameters {
  real<lower=0> T;         // 鶴の数（未知）
  real<lower=0> K;         // 亀の数（未知）
  real<lower=0> sigma_H;   // 観測誤差(観測者の能力): 頭の数
  real<lower=0> sigma_F;   // 観測誤差(観測者の能力): 足の数
  real<lower=2> nu;        // t自由度
}
transformed parameters {
  vector[N] mu_H = rep_vector(T + K, N);        // E[H] = T + K
  vector[N] mu_F = rep_vector(2*T + 4*K, N);    // E[F] = 2T + 4K
}
model {
  (T + K) ~ normal(H_bar, 10);    // 総数は観測された頭の数に近いはず
  sigma_H ~ student_t(3, 0, 10);
  sigma_F ~ student_t(3, 0, 10);
  nu ~ gamma(2, 0.1);             // 平均≈20のロバスト化

  H ~ student_t(nu, mu_H, sigma_H);  // 頭の数 H は mu_H := T + K を中心とした分布に従う
  F ~ student_t(nu, mu_F, sigma_F);  // 足の数 F は mu_F := 2*T + 4*K を中心とした分布に従う
}

データ

冒頭の問題文では "中略" となっていたデータの全貌をお見せしよう。

observed_heads	observed_feet
96	275
96	272
94	274
95	272
97	273
102	266
98	267
100	277
94	268
97	273
95	270
97	272
97	268
96	271
94	265
98	281
95	272
92	275
98	274
96	273

全部で20件ある。

サンプリングと推定

というわけで事後分布を求めて各パラメータを推定した。

	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
T	56.630	54.294	58.854	1.159	6817.0	1.0
K	39.659	38.261	41.033	0.704	6813.0	1.0
sigma_H	2.221	1.423	3.145	0.446	9629.0	1.0
sigma_F	3.836	2.488	5.357	0.756	9835.0	1.0
nu	20.302	2.333	47.689	13.800	9759.0	1.0

鶴の数 T は、平均が 56.630、おおよその範囲が 55\~59 となった。鶴の数 K は、平均が 39.659、おおよその範囲が 38\~41 となった。

「T ≒ 56.630 → 鶴は57羽」のように点推定するだけでなく、確信度に応じて幅を持たせて推定できるのがベイズ推定のうれしいところだ。

また sigma_H, sigma_F の値にも着目してほしい。この値は H, F を観測したときの標準偏差 (ばらつき) を表している。
つまり sigma_H, sigma_F が大きいほど「ばらつきが大きい」=「観測者の能力が低い」と解釈できる。

データを増やして再推定する

追加で20件のデータを取ったとして、推定結果がどう変わるか見てみよう。
追加のデータはこれだ。

observed_heads	observed_feet
95	279
95	273
95	267
94	274
95	271
90	277
95	274
97	277
94	273
95	268
92	270
96	269
94	270
93	277
100	275
95	265
96	270
94	277
97	272
96	272

最初の20件と合わせて合計40件のデータを、先ほどと同じモデルに渡して再度推定してみる。

	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
T	55.071	53.526	56.579	0.776	6995.0	1.0
K	40.513	39.576	41.430	0.474	7003.0	1.0
sigma_H	2.097	1.554	2.723	0.299	9471.0	1.0
sigma_F	3.699	2.810	4.664	0.481	10348.0	1.0
nu	20.076	2.697	45.919	13.147	9186.0	1.0

T, F の平均値はさほど変わっていないが、幅を持たせて推定したときの幅が狭くなっているのに着目してほしい。これはデータを増やしたことで、推定結果により確信が持てたということだ。

	データ20件	データ40件
T の平均	56.630	55.071
T の範囲	54\~59	54\~57
F の平均	39.659	40.513
F の範囲	38\~41	40\~41

答え合わせ

今回使ったデータは私が適当に生成したものなので、実は真の値を知っている。
真の値は T=56, K=40 だ。

	真の値	点推定	区間推定
T	56	55.071	54\~57
F	40	40.513	40\~41

まぁ、いい感じ...？

まとめ

ベイズの考え方を使って、方程式を解くことなく鶴亀算に取り組んでみた。

私からは以上です。

2025-08-26

t分布を使ったロバスト推定

統計 Stan

分析において、得られたデータを分布に当てはめて母数を推定することがよくある。

例えばコイン投げでイカサマコインかどうか見抜きたかったら、とりあえず二項分布に当てはめて「表が出る確率 p」を推定したくなるだろう。

今回は正規分布・t分布への当てはめをやってみる。特に データに外れ値を含む場合にt分布への当てはめが有用 なことを見ていく。
この記事のタイトルにある「ロバスト」とは「t分布を想定しておけば、外れ値を含んだデータに対してもいい感じにやってくれる」みたいな意味だ。

例によって、すべての処理を実行してみられる Google Colab Notebook がここにある。

colab.research.google.com

データ

今回扱うデータを見よう。

[8, 9, 10, 11, 12, 1000]

これだ。めちゃくちゃに単純なデータにしておいた。

大部分は 10 の周りに散らばっている。そしてクソデカ外れ値の 1000 が最後にくっついている。
この 1000 がいかに推定を邪魔するかを見ていく。

モデリング

今回も Stan を使って母数の推定をする。なので Stan コードで「正規分布への当てはめ」と「t分布への当てはめ」を書いておく。

正規分布を使ったモデル

data {
  int<lower=1> N;
  vector[N] X;
}
parameters {
  real mu;              // 平均
  real<lower=0> sigma;  // 標準偏差
}
model {
  X ~ normal(mu, sigma);  // データを正規分布に当てはめる

  mu ~ normal(0, 1000);
  sigma ~ lognormal(0, 0.5);
}

X ~ normal(mu, sigma) としてデータ X が正規分布 (Normal Distribution) に従うことを表現した。
推定する対象は母数 mu, sigma だ。

mu = $\mu$ = 平均
sigma = $\sigma$ = 標準偏差

分布 nomal(mu, sigma) については Stan のリファレンスを見てもらうのが手っ取り早い。

mc-stan.org

t分布を使ったモデル

data {
  int<lower=1> N;
  vector[N] X;
}
parameters {
  real<lower=2, upper=60> nu;  // 自由度
  real mu;                     // 平均
  real<lower=0> sigma;         // 標準偏差
}
model {
  X ~ student_t(nu, mu, sigma);  // データをt分布に当てはめる

  nu ~ gamma(2, 0.2);
  mu ~ normal(0, 1000);
  sigma ~ lognormal(0, 0.5);
}

X ~ student_t(nu, mu, sigma) としてデータ X がt分布 (Student-T Distribution) に従うことを表現した。
正規分布にはなかった母数 nu がある。

nu = $\nu$ = 自由度
mu = $\mu$ = 中心位置
sigma = $\sigma$ = 尺度・スケール

となっている。

student_t(nu, mu, sigma) についても詳細は Stan のリファレンスを見てもらいたい。

mc-stan.org

ちなみに

正規分布の $\mu$ とt分布の $\mu$ は厳密には別物なんだろうけど、「t分布の "中心位置" は正規分布でいう "平均" ですよ」という方便は許されるはずなので同じ記号を当てている。 $\sigma$ についても同様。

まずは外れ値なしで当てはめてみる

まずは外れ値を除いたデータで分布への当てはめをやってみて様子を見よう。

データはこれ。

[8, 9, 10, 11, 12]  # 外れ値 1000 は含めないでおく

正規分布への当てはめ

X = [8, 9, 10, 11, 12]
fit_normal = model_normal.sample(
    data={ "N": len(X), "X": X },
    chains=4, iter_warmup=1000, iter_sampling=10000,
    seed=20250824, adapt_delta=0.9, max_treedepth=12
)

idata_normal = az.from_cmdstanpy(fit_normal)
az.summary(
    idata_normal,
    var_names=["mu", "sigma"],
    hdi_prob=0.95
)[["mean","hdi_2.5%","hdi_97.5%","sd","ess_bulk","r_hat"]]

こんな感じにサンプリングして、推定された母数をサマリーしてみる。

index	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
mu	10.003	8.593	11.421	0.705	18397.0	1.0
sigma	1.502	0.796	2.348	0.426	18311.0	1.0

推定された mu の平均が 10.003 となっている。標本平均 (8+9+10+11+12)/5 → 10 とも整合している。

実際には推定された母数は幅 (分布) を持っている。プロットしたものを見ると mu のおおよその範囲は 8.7~11 と見て取れる。

t分布への当てはめ

同じようにt分布に当てはめて母数を推定する。

X = [8, 9, 10, 11, 12]
fit_student_t = model_student_t.sample(
    data={ "N": len(X), "X": X },
    chains=4, iter_warmup=1000, iter_sampling=10000,
    seed=20250824, adapt_delta=0.9, max_treedepth=12
)

idata_student_t = az.from_cmdstanpy(fit_student_t)
az.summary(
    idata_student_t,
    var_names=["nu", "mu", "sigma"],
    hdi_prob=0.95
)[["mean","hdi_2.5%","hdi_97.5%","sd","ess_bulk","r_hat"]]

index	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
nu	10.727	2.001	24.467	6.974	31269.0	1.0
mu	10.007	8.571	11.519	0.744	27551.0	1.0
sigma	1.406	0.681	2.286	0.437	28410.0	1.0

推定された mu は平均 10.007 で、おおよその範囲は 8.6~11 となっている。

推定された母数で確率密度関数を描く

nu, mu, sigma の値を具体的に推定できたので、正規分布・t分布の式に当てはめて確率密度関数 (PDF) のグラフを描いてみよう。

データ [8, 9, 10, 11, 12] が正規分布・t分布からの実現値と仮定すると、このような形の分布が推定される。
ここまではモデルを正規分布にしてもt分布にしても似たような結果だ。しかしデータに外れ値を含めると違いが見えてくる。

外れ値を含むデータで当てはめ

やっと本題。外れ値を含むデータで当てはめをやってみよう。再度データを見せよう。

[8, 9, 10, 11, 12, 1000]

強烈な外れ値である 1000 が含まれている。

試しに標本平均を計算してみると (8 + 9 + 10 + 11 + 12 + 1000)/6 → 175 となる。単純に標本平均で推定すると直感と異なる結果が得られてしまう。

正規分布への当てはめ

まずは正規分布への当てはめを見よう。

index	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
mu	173.786	27.75	321.778	74.785	23862.0	1.0
sigma	182.94	138.074	235.273	25.487	26742.0	1.0

推定された mu の平均は 173.786、おおよその範囲は 29~311 となった。
外れ値を除いた場合に比べて、推定された mu の範囲が激しく広くなってしまった。

		mean	hdi_2.5%	hdi_97.5%	sd
外れ値除外	mu	10.003	8.593	11.421	0.705
外れ値含む	mu	173.786	27.75	321.778	74.785

t分布への当てはめ

t分布への当てはめでは印象が異なってくる。

index	mean	hdi_2.5%	hdi_97.5%	sd	ess_bulk	r_hat
nu	2.208	2.0	2.62	0.209	25319.0	1.0
mu	10.009	8.16	11.867	0.932	28339.0	1.0
sigma	1.625	0.675	2.842	0.606	28037.0	1.0

推定された mu の平均は 10.009、おおよその範囲は 8.3~12 となった。
なんと外れ値に惑わされず、mu の範囲は 10 前後に留まっている。

		mean	hdi_2.5%	hdi_97.5%	sd
外れ値除外	mu	10.007	8.571	11.519	0.744
外れ値含む	mu	10.009	8.16	11.867	0.932

推定された母数で確率密度関数を描く

再び母数が推定できたので確率密度関数 (PDF) をプロットしてみよう。

外れ値データのせいで正規分布がほぼ平らに潰れてしまっている。

これではよくわからないのでもっと範囲を広げてプロットしてみよう。

めちゃくちゃに潰れてしまっているものの、たしかに x=173 あたりに山があることが分かる。

潰れてしまった正規分布と潰れなかったt分布

推定結果の差がモデルの分布の違いに依るのはあきらかだろう。

正規分布モデルではデータに外れ値 1000 が含まれていることを重く受け止め、「この分布は横に広ーーい形なんだ」という推定になった。
しかしこれは「広い範囲からいろいろな値が得られる」と言っているにすぎず、分析としてはイマイチだ。

一方、t分布では外れ値 1000 の影響は控えめに評価されて、あくまでも「10 周辺の値が得られる」と言っている。
こちらの方が意見がはっきりしていると感じるだろう。

t分布はなぜ外れ値に惑わされない？

ベイズ推定の考え方では、与えられたデータをもとに分布の母数が尤もらしく (最もらしく) なるように調整される。
今回は全ての母数を推定対象としたので、外れ値 1000 もいずれかの母数に影響している。

正規分布の母数は、平均 mu と標準偏差 sigma だ。
外れ値 1000 の影響を反映させようと思ったら、mu と sigma の値をバカでかく調整するしかない。

t分布の場合は自由度 nu, 平均 mu, 標準偏差 sigma を母数に持っている。
自由度 nu が小さくなると分布の裾が厚くなる。自由度 nu が小さいt分布ほど外れ値が出やすいということだ。

Skbkekas - t分布の確率密度関数, CC 表示 3.0, リンクによる

モデルにt分布を使った推定では、外れ値に対して「自由度 nu を小さくする」という調整が可能だ。代わりに mu と sigma の値をそれほど変えずにすむ。
結果として外れ値に惑わされず、mu の推定値は多数派の値のまわりに張り付くことになる。

これがt分布が外れ値に強い理由であり。モデリングでt分布を想定するとロバストになる (外れ値を渡されても惑わされにくくなる) ということだ。

自由度 nu で調整とは

mu と sigma の値を変えずに自由度 nu の値で調整とはどういうことだろうか。t分布の性質を知っておくと理解の助けになる。

nu=1 のときのt分布には「コーシー分布」という特別な名前がついている。そしてコーシー分布には「平均・分散が存在しない」という恐ろしい性質がある。

コーシー分布にも山になっている部分はある (故にその周辺の値が得られやすい) ものの、裾が厚くて外れ値が出やすいために平均が収束しない。
t分布の自由度を nu=1 まで下げるとこのような 性質の悪い分布 になる。

一方、自由度 nu を大きくしていくとt分布は正規分布に近づいていく。nu → ∞ の極限では正規分布に一致する。
正規分布は中心極限定理から自然に導出される分布で、 最も性質のよい分布の一つ と言える。

このようにt分布の自由度 nu は、コーシー分布と正規分布のあいだをグラデーションで繋ぐパラメーターになっている。
nu の値を小さくすることで「外れ値が出やすい分布」を表現できる。

ちなみにt分布の平均が定義されるには 1 < nu でなければいけないし、分散が定義されるには 2 < nu でなければいけない。データ分析する立場からすると 2 < nu であってほしい。

外れ値の有無による nu の変化を見てみる。

というわけでt分布の自由度 nu は「分布の性質の良さ」や「外れ値の出やすさ」を測る目安になる。

nu が、
- 小さい → 性質: 悪い, 外れ値: 多
- 大きい → 性質: 良い, 外れ値: 少

先ほどの推定で得た nu を外れ値の有無によって比較してみよう。

		mean	hdi_2.5%	hdi_97.5%	sd
外れ値除外	nu	10.727	2.001	24.467	6.974
外れ値含む	nu	2.208	2.0	2.62	0.209

除外の方では 10.727 (2.001~24.467) だったのが、含む方では 2.208 (2.0~2.62) まで悪化している。

実際のデータ分析でt分布への当てはめをするときにも、推定された nu の値を見て「正規分布からどれほど離れているのか」を測っていいだろう。

まとめ

モデルにt分布を使うことで外れ値に惑わされにくい推定が可能であることを見た。

私からは以上です。

無駄と文化

実用的ブログ

「決意を新たにする」のはやめよう

「決意を新たにする」という脳内完結型イベント

「辞めないこと」は軽率に再開すること

まとめ

ベイジアン鶴亀算/t分布モデル

ベイズの考えを使って解く

実行環境・ツール

モデル定義

Stan コード

データ

サンプリングと推定

データを増やして再推定する

答え合わせ

まとめ

t分布を使ったロバスト推定

データ

モデリング

正規分布を使ったモデル

t分布を使ったモデル

ちなみに

まずは外れ値なしで当てはめてみる

正規分布への当てはめ

t分布への当てはめ

推定された母数で確率密度関数を描く

外れ値を含むデータで当てはめ

正規分布への当てはめ

t分布への当てはめ

推定された母数で確率密度関数を描く

潰れてしまった正規分布と潰れなかったt分布

t分布はなぜ外れ値に惑わされない？

自由度 nu で調整とは

外れ値の有無による nu の変化を見てみる。

まとめ

observed_heads	observed_feet
96	275
96	272
94	274
95	272
97	273
102	266
98	267
100	277
94	268
97	273
95	270
97	272
97	268
96	271
94	265
98	281
95	272
92	275
98	274
96	273

observed_heads	observed_feet
95	279
95	273
95	267
94	274
95	271
90	277
95	274
97	277
94	273
95	268
92	270
96	269
94	270
93	277
100	275
95	265
96	270
94	277
97	272
96	272

observed_heads	observed_feet
96	275
96	272
94	274
95	272
97	273
102	266
98	267
100	277
94	268
97	273
95	270
97	272
97	268
96	271
94	265
98	281
95	272
92	275
98	274
96	273

observed_heads	observed_feet
95	279
95	273
95	267
94	274
95	271
90	277
95	274
97	277
94	273
95	268
92	270
96	269
94	270
93	277
100	275
95	265
96	270
94	277
97	272
96	272

observed_heads	observed_feet
96	275
96	272
94	274
95	272
97	273
102	266
98	267
100	277
94	268
97	273
95	270
97	272
97	268
96	271
94	265
98	281
95	272
92	275
98	274
96	273

observed_heads	observed_feet
95	279
95	273
95	267
94	274
95	271
90	277
95	274
97	277
94	273
95	268
92	270
96	269
94	270
93	277
100	275
95	265
96	270
94	277
97	272
96	272