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分布を使うことで外れ値に惑わされにくい推定が可能であることを見た。

私からは以上です。

2025-08-25

【SQL】期間を指定して日付を列挙する

BigQuery SQL

モチベーション

SQL を書いていると期間をもとに日付を列挙したくなることが稀によくある。書き方を見てみよう。

免責事項

本記事は BigQuery を前提にしています。紹介している generate_date_array 関数は MySQL/PostgreSQL には存在しない BigQuery 独自の機能です。
ちなみに PostgreSQL では generate_series 関数を使って同様のことができるらしいです。

`unnest(generate_date_array(...))` を使おう

手元に START_DATE, END_DATE の値があるとして、

select d
from unnest(generate_date_array(START_DATE, END_DATE)) as d

こう書くとよい。

試そう。

with
_dates as (
    select d
    from unnest(generate_date_array('2025-08-10', '2025-08-13'))
        as d
)

select *
from _dates
order by d asc

2025-08-10 と 2025-08-13 を指定して、

	d
1	2025-08-10
2	2025-08-11
3	2025-08-12
4	2025-08-13

が得られた

応用: 特定の曜日だけ列挙する

WHERE 句と組み合わせると条件を満たす日付だけを列挙できる。
例えば 2025年8月中の火曜, 金曜がほしいなら、

with
_dates as (
    select d
    from unnest(generate_date_array('2025-08-01', '2025-08-31'))
        as d
    where format_date('%a', d) in ('Tue', 'Fri')
)

select d, format_date('%a', d) as day_of_week
from _dates
order by d asc

こうだ。

WHERE 句に where format_date('%a', d) in ('Tue', 'Fri') を指定しているので火曜, 金曜だけが抽出される。

2025-08-01 と 2025-08-31 を指定して、

	d	day_of_week
1	2025-08-01	Fri
2	2025-08-05	Tue
3	2025-08-08	Fri
4	2025-08-12	Tue
5	2025-08-15	Fri
6	2025-08-19	Tue
7	2025-08-22	Fri
8	2025-08-26	Tue
9	2025-08-29	Fri

が得られた。

起こっていることを理解しよう

GENERATE_DATE_ARRAY 関数は開始日と終了日を渡すと ARRAY 型の値を返してくれる。

select generate_date_array('2025-08-10', '2025-08-13');
/*----------------------------------------------------------*
 | f0_                                                      |
 +----------------------------------------------------------+
 | ['2025-08-10', '2025-08-11', '2025-08-12', '2025-08-13'] |
 *----------------------------------------------------------*/

UNNEST 演算子は ARRAY をフラット化して行の集合に変換してくれる。

select *
from unnest(generate_date_array('2025-08-10', '2025-08-13'));
/*--------------*
 | f0_          |
 +--------------+
 | '2025-08-10' |
 | '2025-08-11' |
 | '2025-08-12' |
 | '2025-08-13' |
 *--------------*/

N行×1列のテーブル相当の結果が返るので、UNNEST はサブクエリが書ける場所でしか使えない。

unnest(...) が関数っぽく見えてしまうので、結果として1列のテーブルが返るのは相当混乱する。
実際のクエリでは CTE を使って、

with

_dates as (
    select d
    from unnest(generate_date_array('2025-08-10', '2025-08-13')) as d;
)

...

と書いてしまって。以降は _dates というテーブルとして扱うのが混乱が少ないと思う。

2025-08-17

StanによるMCMCサンプリング/素振り/ベータ分布

統計 Stan

MCMC サンプリングの練習をしよう。
Stan を使う。実行環境は Google Colab 、よってローカル環境でのセットアップは一切不要。

colab.research.google.com

ここにノートブックがある。

モデル定義

シンプルにただのベータ分布からサンプリングするだけをやってみよう。

data {
    real<lower=0> alpha;
    real<lower=0> beta;
}
parameters {
    real<lower=0, upper=1> theta;
}
model {
    theta ~ beta(alpha, beta);
}

やってることは $\theta \sim beta(\alpha, \beta)$ だけ。この $\theta$ をサンプリングして評価すればいい。

サンプリング実行

モデルをコンパイルしてサンプリングを実行しよう。

beta_model = CmdStanModel(stan_file=str(model_path))

fit = beta_model.sample(
    data={"alpha": 2.0, "beta": 5.0},
    chains=8, parallel_chains=8,
    iter_warmup=1000, iter_sampling=10000,
    seed=20250815
)

$\alpha = 2, \beta = 5$ として $beta(\alpha = 2, \beta = 5)$ からサンプリングする。
8並列で10,000個ずつサンプルを取る。最初の1,000個はバーンインとして捨てる。最終的に80,000個のサンプルが得られる。

サンプリング結果

得られたサンプルでヒストグラムを描いてみよう。

いい感じにフィットして見える。

サンプルから統計量を推測する

得られたサンプルから標本平均, 標本標準偏差を計算すると、真の平均や真の標準偏差を推測できる。

summ = fit.summary()
print(summ.loc[["theta"], ["Mean", "StdDev", "R_hat", "ESS_bulk", "ESS_tail"]])

結果は、

	Mean (平均)	StdDev (標準偏差)	R_hat	ESS_bulk	ESS_tail
theta	0.284623	0.159601	1.00022	26463.6	25729.0

となった。

ベータ分布の平均と標準偏差の理論値は、

$\mu = \frac{\alpha}{\alpha + \beta}$
$\sigma = \frac{\alpha \beta}{(\alpha + \beta)^{2}(\alpha + \beta + 1)}$

なのでここに $\alpha = 2, \beta = 5$ を代入して比較すると、

	平均	標準偏差
理論値	0.285714	0.159718
標本からの推定	0.284623	0.159601

まぁまぁ良さそう。

まとめ

Google Colab は便利。Stan の文法に慣れていきたい。

私からは以上です。

無駄と文化

実用的ブログ

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

データ

モデリング

正規分布を使ったモデル

t分布を使ったモデル

ちなみに

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

正規分布への当てはめ

t分布への当てはめ

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

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

正規分布への当てはめ

t分布への当てはめ

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

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

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

自由度 nu で調整とは

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

まとめ

【SQL】期間を指定して日付を列挙する

モチベーション

免責事項

`unnest(generate_date_array(...))` を使おう

応用: 特定の曜日だけ列挙する

起こっていることを理解しよう

StanによるMCMCサンプリング/素振り/ベータ分布

モデル定義

サンプリング実行

サンプリング結果

サンプルから統計量を推測する

まとめ

データ

モデリング

正規分布を使ったモデル

t分布を使ったモデル

ちなみに

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

正規分布への当てはめ

t分布への当てはめ

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

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

正規分布への当てはめ

t分布への当てはめ

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

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

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

自由度 nu で調整とは

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

まとめ

モチベーション

免責事項

unnest(generate_date_array(...)) を使おう

応用: 特定の曜日だけ列挙する

起こっていることを理解しよう

モデル定義

サンプリング実行

サンプリング結果

サンプルから統計量を推測する

まとめ

`unnest(generate_date_array(...))` を使おう