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

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

例えばコイン投げでイカサマコインかどうか見抜きたかったら、とりあえず二項分布に当てはめて「表が出る確率 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分布を使ったロバスト推定

データ

モデリング

正規分布を使ったモデル

t分布を使ったモデル

ちなみに

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

正規分布への当てはめ

t分布への当てはめ

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

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

正規分布への当てはめ

t分布への当てはめ

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

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

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

自由度 nu で調整とは

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

まとめ