手元に300個の測定値があって、「これはガンマ分布だろう」と当たりはついている。では形状パラメータと尺度パラメータをいくらにすればデータにフィットするのでしょうか。尤度関数を立てて対数を取り、偏微分してゼロと置いて……と最尤推定の手続きを踏むこともできます。しかしガンマ分布の最尤方程式にはディガンマ関数が現れ、手では解けず数値解法が要ります。
ところが、もっと素朴で速い方法があります。「データの平均と分散を計算して、それがガンマ分布の理論的な平均・分散と一致するように、パラメータを逆算する」だけ。これが モーメント法（method of moments、積率推定） です。19世紀末にカール・ピアソンが体系化した、推定法の中でも最も古く、そして最も直感的な手法です。
モーメント法は今でも現役です。たとえば、保険数理や信頼性工学で観測データに分布をフィットする「分布フィッティング」の第一手として。最尤推定の反復計算を始める前の 初期値 を一瞬で得る道具として。そして計量経済学では、モーメント法を一般化した 一般化モーメント法（GMM） が、内生性をもつ回帰モデルの推定における標準ツールになっています。簡単なのに応用が広い——それがモーメント法です。
本記事の内容
- 標本モーメントと理論モーメントを「等置する」という直感
- 推定量を求める一般的な手順（k個のパラメータにはk個のモーメント）
- 正規・ガンマ・ベータ・一様分布での具体的な導出
- 最尤推定（MLE）との比較 — 簡単だが効率は劣る、しかし初期値として有用
- 一致性と漸近正規性（大数の法則とデルタ法）
- 計量経済学への発展としての一般化モーメント法（GMM）
- Pythonによる複数分布フィットの実装と考察
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
モーメントを「合わせる」という直感
体重計を思い浮かべてください。あなたは体重計の目盛りがどこかでズレている（バイアスがある）かもしれないと疑っている。確かめるには、すでに重さがわかっている標準分銅を載せて、表示が分銅の本当の重さと一致するように目盛りを調整しますよね。「測った値」を「あるべき値」に合わせて、機器のパラメータを決める——これがモーメント法の発想そのものです。
統計の文脈では、「測った値」がデータから計算する 標本モーメント、「あるべき値」が分布の理論から決まる 理論モーメント です。$k$ 次のモーメントとは、ざっくり言えば「$k$ 乗の平均」のこと。1次モーメントは平均、2次モーメント（中心化したもの）は分散に対応します。
データ $x_1, x_2, \dots, x_n$ から計算する $k$ 次の標本モーメントは
$$ \hat{m}_k = \frac{1}{n}\sum_{i=1}^{n} x_i^{\,k} $$
です。一方、パラメータ $\theta$ をもつ分布の $k$ 次の理論モーメントは、期待値
$$ \mu_k(\theta) = E[X^k] $$
で与えられ、$\theta$ の関数になります。モーメント法のアイデアは、この2つを 等置（イコールで結ぶ） することです。
$$ \hat{m}_k = \mu_k(\theta) $$
下の図が、この発想の全体像です。
左の「データから計算する量」と、右の「パラメータの関数として表される量」を、真ん中で等号で結ぶ。あとはその方程式を $\theta$ について解けば推定量が得られます。求めたいパラメータが2個なら、1次と2次の2本の方程式を連立させればよい。なぜこれでうまくいくのかというと、大数の法則 により標本モーメント $\hat{m}_k$ はサンプル数が増えれば理論モーメント $\mu_k(\theta_0)$（$\theta_0$ は真のパラメータ）に近づくからです。つまり「測った値 ≈ あるべき値」が漸近的に保証されている。
この素朴な等置がなぜ妥当なのかが見えたところで、実際にパラメータを求める手順を整理しましょう。
モーメント法の手順
手順は驚くほど機械的です。パラメータが $p$ 個あるとき、次の3ステップで推定量が求まります。
ステップ1: 理論モーメントをパラメータの式で書く。 推定したいパラメータが $\theta = (\theta_1, \dots, \theta_p)$ の $p$ 個なら、1次から $p$ 次までの理論モーメント $\mu_1(\theta), \dots, \mu_p(\theta)$ を、$\theta$ の関数として書き下します。多くの教科書的分布では、これらは公式として知られています（正規なら $\mu_1 = \mu$、$\mu_2 = \mu^2 + \sigma^2$ など）。
ステップ2: 標本モーメントと等置する。 各次数で、標本モーメント $\hat{m}_k$ と理論モーメント $\mu_k(\theta)$ をイコールで結び、$p$ 本の連立方程式を作ります。
$$ \begin{aligned} \mu_1(\theta_1, \dots, \theta_p) &= \hat{m}_1 \\ \mu_2(\theta_1, \dots, \theta_p) &= \hat{m}_2 \\ &\ \ \vdots \\ \mu_p(\theta_1, \dots, \theta_p) &= \hat{m}_p \end{aligned} $$
ステップ3: 連立方程式を $\theta$ について解く。 解 $\hat{\theta} = (\hat{\theta}_1, \dots, \hat{\theta}_p)$ がモーメント法推定量（MoM推定量）です。
なお、平均まわりの量（分散など）を扱うときは、生のモーメント $\hat{m}_k$ の代わりに 中心モーメント を使うと式が簡単になることがよくあります。2次の中心モーメント（標本分散）は
$$ \hat{\mu}_2^{(c)} = \frac{1}{n}\sum_{i=1}^n (x_i – \bar{x})^2 = \hat{m}_2 – \hat{m}_1^2 $$
です。ここで割る数は $n-1$ ではなく $n$ であることに注意してください。モーメント法は「標本モーメントをそのまま使う」のが定義なので、不偏分散（$n-1$ で割る）ではなく 標本分散（$n$ で割る） を用います。本記事の数式・コードはこの $1/n$ の流儀で統一します。
たったこれだけです。最尤推定のように尤度関数を微分する必要はありません。では、この手順を実際の分布に当てはめてみましょう。
具体例1: 正規分布
正規分布 $N(\mu, \sigma^2)$ のパラメータは $\mu$ と $\sigma^2$ の2個なので、1次と2次の2本のモーメントを使います。
理論モーメントは、正規分布の定義から次のとおりです。1次モーメント（平均）はそのまま $\mu$、2次の生モーメントは「$E[X^2] = E[X]^2 + V[X]$」より $\mu^2 + \sigma^2$ です。
$$ \mu_1 = E[X] = \mu, \qquad \mu_2 = E[X^2] = \mu^2 + \sigma^2 $$
これを標本モーメントと等置します。
$$ \hat{m}_1 = \mu, \qquad \hat{m}_2 = \mu^2 + \sigma^2 $$
1本目から直ちに $\hat{\mu} = \hat{m}_1 = \bar{x}$ が得られます。これを2本目に代入して $\sigma^2$ について解きましょう。
$$ \hat{\sigma}^2 = \hat{m}_2 – \hat{\mu}^2 = \hat{m}_2 – \hat{m}_1^2 = \frac{1}{n}\sum_i x_i^2 – \bar{x}^2 $$
最後の式は、まさに先ほどの標本分散 $\frac{1}{n}\sum_i (x_i – \bar{x})^2$ に他なりません。つまり正規分布のMoM推定量は「平均は標本平均、分散は標本分散」という、誰もが直感的にやる推定そのものです。
実際にデータでやってみたのが次の図です。
この図から2つのことが読み取れます。第一に、オレンジのMoM推定PDFが青破線の真の分布とほぼ重なっており、平均・分散を合わせるだけで分布全体がよく再現できています。第二に、推定された $\hat{\mu} = 2.96$、$\hat{\sigma}^2 = 1.96$ は真値 $\mu = 3.0$、$\sigma^2 = 1.96$ にきわめて近い。正規分布の場合、実はMoM推定量は（分散の分母の $n$ vs $n-1$ の違いを除けば）最尤推定量と一致します。この特別な一致は、後で見るように一般には成り立ちません。
正規分布は最も簡単なケースでした。次は、最尤推定が手では解けないガンマ分布で、モーメント法の威力を見ましょう。
具体例2: ガンマ分布
ガンマ分布は形状パラメータ $k$ と尺度パラメータ $\theta$ をもち、待ち時間や正の量のモデルに広く使われます。この分布の最尤推定はディガンマ関数を含む方程式になり、数値解法が必要です。ところがモーメント法なら、紙とペンだけで閉じた式が得られます。
ガンマ分布 $\mathrm{Gamma}(k, \theta)$ の理論的な平均と分散は
$$ E[X] = k\theta, \qquad V[X] = k\theta^2 $$
です。標本平均 $\bar{x}$ と標本分散 $\hat{\sigma}^2$ にそれぞれ等置します。
$$ k\theta = \bar{x}, \qquad k\theta^2 = \hat{\sigma}^2 $$
ここで $\theta$ を消すために、2本目を1本目で割ります。すると左辺は $\frac{k\theta^2}{k\theta} = \theta$ となり、尺度パラメータがすぐ求まります。
$$ \hat{\theta} = \frac{k\theta^2}{k\theta} = \frac{\hat{\sigma}^2}{\bar{x}} $$
この $\hat{\theta}$ を1本目の $k\theta = \bar{x}$ に戻して $k$ について解くと、
$$ \hat{k} = \frac{\bar{x}}{\hat{\theta}} = \frac{\bar{x}}{\hat{\sigma}^2 / \bar{x}} = \frac{\bar{x}^2}{\hat{\sigma}^2} $$
が得られます。形状パラメータは「平均の2乗 ÷ 分散」、尺度パラメータは「分散 ÷ 平均」。どちらも電卓があれば一瞬で計算できる、見事に閉じた式です。
実データへのフィットが次の図です。
推定値は $\hat{k} = 2.33$、$\hat{\theta} = 1.70$ で、真値 $k = 2.5$、$\theta = 1.6$ に近い値が、反復計算なしで得られています。オレンジのMoM推定曲線はヒストグラムの山の位置と裾の広がりをよく捉えています。ピアソンがモーメント法を編み出した最大の動機が、まさにこうした「最尤推定が解析的に解けない分布」を手早く扱うことだったのです。
ガンマ分布では平均・分散を割り算するだけでした。次のベータ分布では、もう少し代数の操作が要りますが、考え方は同じです。
具体例3: ベータ分布
ベータ分布 $\mathrm{Beta}(\alpha, \beta)$ は区間 $[0, 1]$ 上の分布で、確率や割合のモデルに使われます。理論的な平均と分散は次のとおりです。
$$ E[X] = \frac{\alpha}{\alpha + \beta}, \qquad V[X] = \frac{\alpha\beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)} $$
標本平均を $\bar{x} = m$、標本分散を $v$ とおいて等置します。まず平均の式 $m = \frac{\alpha}{\alpha+\beta}$ から $1 – m = \frac{\beta}{\alpha+\beta}$ なので、この2つを掛けると
$$ m(1 – m) = \frac{\alpha}{\alpha+\beta} \cdot \frac{\beta}{\alpha+\beta} = \frac{\alpha\beta}{(\alpha+\beta)^2} $$
という関係が出ます。これは分散の式の分子・分母にそっくりです。分散の式を、いま得た関係を使って書き直してみましょう。分散の式の分子 $\alpha\beta$ と分母の $(\alpha+\beta)^2$ をまとめると、
$$ v = \frac{\alpha\beta}{(\alpha+\beta)^2 (\alpha+\beta+1)} = \frac{m(1-m)}{\alpha + \beta + 1} $$
となります。$\frac{\alpha\beta}{(\alpha+\beta)^2}$ を $m(1-m)$ で置き換えただけです。この式を $\alpha + \beta$ について解くと、
$$ \alpha + \beta = \frac{m(1-m)}{v} – 1 $$
が得られます。右辺はすべてデータから計算できる量です。この「合計」を $S = \frac{m(1-m)}{v} – 1$ とおきましょう。あとは $\alpha = m(\alpha + \beta) = mS$、$\beta = (1-m)(\alpha+\beta) = (1-m)S$ と、平均の式から各パラメータを取り出せます。
$$ \hat{\alpha} = m\left(\frac{m(1-m)}{v} – 1\right), \qquad \hat{\beta} = (1-m)\left(\frac{m(1-m)}{v} – 1\right) $$
これでベータ分布のMoM推定量が閉じた形で求まりました。実装結果が次の図です。
推定値 $\hat{\alpha} = 1.97$、$\hat{\beta} = 4.95$ は真値 $\alpha = 2$、$\beta = 5$ にぴったり近く、左に偏った非対称な形がよく再現されています。ベータ分布の最尤推定もディガンマ関数を含んで数値解法が必要なので、ここでもモーメント法の「閉じた式で一発」という利点が活きています。
ここまでの3例は「中心極限的な滑らかな分布」でした。次の一様分布は、モーメント法と最尤推定の性格の違いがはっきり出る、教訓的な例です。
具体例4: 一様分布
連続一様分布 $U(a, b)$ は区間 $[a, b]$ 上で一定の密度をもちます。パラメータは下端 $a$ と上端 $b$ の2個。理論的な平均と分散は次のとおりです。
$$ E[X] = \frac{a + b}{2}, \qquad V[X] = \frac{(b – a)^2}{12} $$
平均の式から $a + b = 2m$（$m$ は標本平均）。分散の式から $(b – a)^2 = 12 v$ なので $b – a = 2\sqrt{3}\,\sqrt{v} = 2\sqrt{3}\,s$（$s$ は標本標準偏差）が出ます。$\sqrt{12} = 2\sqrt{3}$ を使いました。和と差が分かれば各々が求まるので、
$$ \hat{a} = m – \sqrt{3}\,s, \qquad \hat{b} = m + \sqrt{3}\,s $$
がMoM推定量です。「標本平均を中心に、標準偏差の $\sqrt{3}$ 倍だけ両側に広げる」というわけです。
一方、一様分布の最尤推定はまったく別物で、$\hat{a} = \min_i x_i$、$\hat{b} = \max_i x_i$、つまり データの最小値と最大値 そのものになります。この2つを比べたのが次の図です。
この図はモーメント法の 弱点 を露わにします。MoM推定の端点（オレンジ実線）は、観測されたデータの範囲を はみ出す ことがあります。図でも左端のMoM推定がデータの最小値より外側に出ています。逆に、はみ出さない場合は端点が内側に来てしまい、「観測されたのに区間外」という矛盾が起き得ます。MLEの端点（緑点線）はデータの最小・最大なので、少なくともすべての観測点を含みます。一様分布のように 台（サポート）の境界がパラメータになっている分布 では、モーメント法は素直に機能しないことを覚えておきましょう。
このように、モーメント法とMLEは結果が一致するとは限りません。両者の優劣を、推定量の「ばらつき」という観点から正面から比較しましょう。
最尤推定との比較
モーメント法と最尤推定（MLE）は、推定の哲学が違います。モーメント法は「いくつかの要約統計量（平均・分散）だけを合わせる」のに対し、MLEは「全データの尤度を最大化する」。MLEはデータの情報を最大限に使うため、一般に 効率的（同じサンプル数でばらつきが小さい）です。漸近的には、MLEは達成可能な最小分散（クラメール・ラオの下界）を達成します。モーメント法はそこまでの保証がありません。
実際にどれくらい差が出るのか、ガンマ分布の形状パラメータ $k$ で見てみましょう。同じサンプル数 $n = 30$ のデータを4000回生成し、それぞれでMoM推定値とMLE推定値を計算して、推定値のヒストグラムを描いたのが次の図です。
この図から、推定法の本質的な違いがはっきり見えます。両者とも分布の中心は真値 $k = 3$ あたりにあり、平均的には正しい値を当てています。しかし ばらつき（分散） が違う。MoM推定（オレンジ）の分散は 1.028、MLE推定（緑）の分散は 0.778 で、MLEの方が3割ほどばらつきが小さい。これが「MLEの方が効率的」ということの具体的な意味です。同じデータからより信頼できる推定値を引き出せる。
では、なぜ効率で劣るモーメント法をわざわざ使うのでしょうか。理由は3つあります。
第一に、計算が簡単で速い。 閉じた式があれば反復計算が不要です。大量の分布を次々にフィットする場面で効きます。
第二に、MLEの初期値として優秀。 ガンマやベータの最尤推定は数値最適化（ニュートン法など）で解きますが、最適化には初期値が要ります。デタラメな初期値だと収束に時間がかかったり、局所解に捕まったりします。モーメント法推定量は真値の近くにあることが保証されている（一致性、後述）ので、最良の初期値 になります。実際、scipy.stats の fit も内部でモーメント法的な初期値を使っています。
第三に、尤度が書けない／複雑すぎる場合でも使える。 モーメント条件さえ立てられれば推定できる。この性質を極限まで一般化したのが、後で扱うGMMです。
「平均的には正しい」と書きましたが、サンプル数を増やせば本当に真値に収束するのでしょうか。次はその理論的保証、一致性を見ます。
一致性と漸近正規性
モーメント法推定量が信頼できる根拠は、2つの定理に支えられています。大数の法則 と 中心極限定理（およびデルタ法）です。
一致性（consistency） とは、サンプル数 $n$ を無限大にすると推定量が真値に確率収束する性質です。なぜモーメント法推定量が一致性をもつのか。出発点は大数の法則です。標本モーメントは、$n \to \infty$ で理論モーメントに確率収束します。
$$ \hat{m}_k = \frac{1}{n}\sum_{i=1}^n x_i^k \ \xrightarrow{\ p\ }\ E[X^k] = \mu_k(\theta_0) $$
ここで $\theta_0$ は真のパラメータです。モーメント法推定量は、理論モーメントを標本モーメントの関数として逆に解いたもの $\hat{\theta} = g(\hat{m}_1, \dots, \hat{m}_p)$ です。この関数 $g$（モーメントからパラメータを復元する写像）が連続なら、連続写像定理により
$$ \hat{\theta} = g(\hat{m}_1, \dots, \hat{m}_p) \ \xrightarrow{\ p\ }\ g(\mu_1(\theta_0), \dots, \mu_p(\theta_0)) = \theta_0 $$
が成り立ちます。「標本モーメントが真のモーメントに近づく」+「モーメントからパラメータへの写像が連続」⇒「推定量が真値に近づく」という三段論法です。
この収束の様子を数値実験で見たのが次の図です。
横軸は対数スケールのサンプル数 $n$、縦軸はガンマ形状パラメータ $k$ の推定値です。$n$ が小さいうちは推定値が大きくばらつき（薄いオレンジの帯が広い）、平均値も真値からずれがちです。しかし $n$ を増やすにつれて、平均値（濃いオレンジ線）が真値 $k = 2.5$（黒破線）に収束し、ばらつきの帯がきれいに縮んでいきます。これが一致性の視覚的な現れです。
さらに、適切な条件のもとでモーメント法推定量は 漸近正規性 をもちます。すなわち $\sqrt{n}(\hat{\theta} – \theta_0)$ が漸近的に正規分布に従う。これは中心極限定理（標本モーメントの漸近正規性）に、写像 $g$ の線形近似を施す デルタ法 を適用すれば導けます。漸近正規性のおかげで、モーメント法推定量に対しても標準誤差や信頼区間を構成できます。ただし、その漸近分散はMLEのものより大きい（効率で劣る）のが一般的でした。
ここまでは「パラメータ数とモーメント条件の数がぴったり等しい」場合でした。では、モーメント条件の方が多い、過剰な情報があるときはどうするか。それが計量経済学を変えたGMMの問いです。
一般化モーメント法（GMM）
これまでのモーメント法は、パラメータ $p$ 個に対してちょうど $p$ 本のモーメント条件を立て、連立方程式を解いてきました。これを「ちょうど識別（just-identified）」と呼びます。ところが現実には、パラメータ数より 多くの モーメント条件が手に入ることがあります。たとえば計量経済学では、ある仮定（操作変数の外生性など）から「この量の期待値はゼロのはずだ」という条件が、パラメータ数より多く湧き出てくる。
モーメント条件を一般的に書くと、ある関数 $g$ について
$$ E[\,\bm{g}(X, \theta_0)\,] = \bm{0} $$
という形になります。ここで $\bm{g}$ は $L$ 本の条件をまとめたベクトル、$\theta$ は $p$ 個のパラメータです。ふつうのモーメント法は $g_k(X, \theta) = X^k – \mu_k(\theta)$ とおいた特別な場合に当たります。
問題は $L > p$ のとき、つまり 過剰識別（over-identified） のときです。条件が多すぎて、すべてを同時にぴったりゼロにする $\theta$ は一般に存在しません。下の図がこの状況です。
ラース・ハンセンが1982年に提案した 一般化モーメント法（Generalized Method of Moments、GMM） の解決策はシンプルです。「全部をゼロにできないなら、全体としてゼロに最も近づける $\theta$ を選ぶ」。具体的には、標本でのモーメント条件の平均
$$ \bar{\bm{g}}(\theta) = \frac{1}{n}\sum_{i=1}^n \bm{g}(x_i, \theta) $$
の大きさを、重み行列 $\bm{W}$ を使った二次形式で測り、それを最小化します。
$$ \hat{\theta}_{\mathrm{GMM}} = \arg\min_{\theta}\ \bar{\bm{g}}(\theta)^{\top}\, \bm{W}\, \bar{\bm{g}}(\theta) $$
$\bm{W}$ は各モーメント条件の重み付けを決める正定値行列です。「精度の高い（分散の小さい）条件を重視する」のが効率的な選び方で、最適な $\bm{W}$ はモーメント条件の共分散行列の逆行列だと知られています。$L = p$ のちょうど識別なら、$\bar{\bm{g}}(\theta) = \bm{0}$ を完全に満たす解が存在するので、$\bm{W}$ によらず通常のモーメント法に一致します。
GMMが計量経済学で絶大な影響をもったのは、尤度関数を書く必要がない からです。誤差項の分布を仮定せず、「期待値に関する条件」だけからパラメータを推定できる。操作変数法、動学パネルデータ、資産価格モデルなど、応用は枚挙にいとまがありません。ハンセンはこの業績で2013年にノーベル経済学賞を受賞しました。
過剰識別のもう一つの利点は モデル検定 です。条件が余っているので、推定後に「余った条件が本当にゼロ近くか」を調べれば、モデルの妥当性を検定できます（ハンセンのJ検定）。情報を捨てずに検証に回せるわけです。
GMMの抽象的な定式化を、もう少し具体的に「条件がゼロになる点を探す」という幾何イメージで掴んでおきましょう。
モーメント条件の幾何イメージ
モーメント法もGMMも、結局は「モーメント条件 $g(\theta) = 0$ を満たす $\theta$ を探す」作業です。これを1次元で可視化してみます。
指数分布 $\mathrm{Exp}(\lambda)$ を例にとります。レートパラメータ $\lambda$ の理論モーメントは $E[X] = 1/\lambda$、$E[X^2] = 2/\lambda^2$ です。それぞれから2つのモーメント条件が作れます。
$$ g_1(\lambda) = \hat{m}_1 – \frac{1}{\lambda}, \qquad g_2(\lambda) = \hat{m}_2 – \frac{2}{\lambda^2} $$
それぞれの条件を $\lambda$ の関数として描き、横軸（ゼロ）を切る点を見たのが次の図です。
$g_1(\lambda) = 0$ の解（青い点）と $g_2(\lambda) = 0$ の解（オレンジの点）は、どちらも真値 $\lambda = 1.5$（黒破線）の近くにありますが、わずかにずれています。これがまさに過剰識別の状況です。パラメータは1個なのに条件が2本あるので、両方を同時にゼロにする $\lambda$ は（有限標本では）存在しない。1次モーメントだけ使えば青い点、2次モーメントだけ使えばオレンジの点が得られ、どちらを採用するかで答えが変わる。
GMMは「2つの条件のどちらか一方を選ぶ」のではなく、「両方の条件をできるだけ小さくする中間点を、精度に応じた重みで選ぶ」という立場をとります。この図を見ると、複数のモーメント条件をどう統合するかという問題が自然に立ち上がってくるのが分かります。GMMはその統合のための一般原理なのです。
理論を一通り見たので、複数の分布を一気にフィットするコードでモーメント法を体感しましょう。
Pythonでの実装
ここまでの分布をモーメント法でフィットする関数を実装します。まずは正規・ガンマ・ベータ・一様の4分布それぞれのMoM推定量を、これまで導いた閉じた式そのままに書き下します。
import numpy as np
from scipy import stats
def mom_normal(x):
"""正規分布: μ=平均, σ²=標本分散(1/n)"""
m1 = x.mean()
var = x.var() # numpyのvarはデフォルトで1/n（ddof=0）
return {"mu": m1, "sigma2": var}
def mom_gamma(x):
"""ガンマ分布: k=平均²/分散, θ=分散/平均"""
m1 = x.mean()
var = x.var()
theta = var / m1
k = m1 / theta # = m1² / var
return {"k": k, "theta": theta}
def mom_beta(x):
"""ベータ分布: α=mS, β=(1-m)S, S=m(1-m)/v - 1"""
m = x.mean()
v = x.var()
S = m * (1 - m) / v - 1
return {"alpha": m * S, "beta": (1 - m) * S}
def mom_uniform(x):
"""一様分布: a=平均-√3·std, b=平均+√3·std"""
m = x.mean()
s = x.std() # 1/n の標準偏差
half = np.sqrt(3) * s
return {"a": m - half, "b": m + half}
4つとも、平均（mean）と分散（var）を計算して代入するだけの、ごく短い関数です。最尤推定のような反復最適化はどこにもありません。これがモーメント法の手軽さです。実際に乱数データで動かして、推定値が真値に近いか確認しましょう。
rng = np.random.default_rng(1)
# 正規分布
x = rng.normal(3.0, 1.4, 400)
print("正規 真値: mu=3.0, sigma2=1.96")
print("正規 MoM:", {k: round(v, 3) for k, v in mom_normal(x).items()})
# ガンマ分布
x = rng.gamma(2.5, 1.6, 500)
print("ガンマ 真値: k=2.5, theta=1.6")
print("ガンマ MoM:", {k: round(v, 3) for k, v in mom_gamma(x).items()})
# ベータ分布
x = rng.beta(2.0, 5.0, 500)
print("ベータ 真値: alpha=2.0, beta=5.0")
print("ベータ MoM:", {k: round(v, 3) for k, v in mom_beta(x).items()})
実行すると、ガンマは k≈2.33, theta≈1.70、ベータは alpha≈1.97, beta≈4.95 のように、いずれも真値の近くの値が一瞬で返ります。サンプル数が数百でこの精度なら、分布フィットの第一手として十分実用的です。乱数のseedにより数値は多少前後しますが、真値の周りに集まる傾向は変わりません。
次に、モーメント法がMLEの 初期値 としていかに優秀かを、ガンマ分布の最尤推定で確かめます。
from scipy.optimize import minimize
from scipy.special import digamma, gammaln
def neg_loglik_gamma(params, x):
"""ガンマ分布の負の対数尤度（形状k, 尺度theta）"""
k, theta = params
if k <= 0 or theta <= 0:
return np.inf
n = len(x)
ll = (k - 1) * np.log(x).sum() - x.sum() / theta \
- n * k * np.log(theta) - n * gammaln(k)
return -ll
rng = np.random.default_rng(2)
x = rng.gamma(2.5, 1.6, 500)
# モーメント法で初期値を作る
init = mom_gamma(x)
x0 = [init["k"], init["theta"]]
# その初期値からMLEを最適化
res = minimize(neg_loglik_gamma, x0, args=(x,), method="Nelder-Mead")
print("MoM初期値:", [round(v, 3) for v in x0])
print("MLE収束値:", [round(v, 3) for v in res.x])
print("最適化の反復回数:", res.nit)
出力を見ると、MoMの初期値 [2.33, 1.70] 付近からMLEがすぐに収束し、[2.4, 1.66] のような最尤推定値に少ない反復回数でたどり着きます。デタラメな初期値（たとえば [100, 0.01]）から始めると反復回数が増えたり収束が不安定になったりしますが、モーメント法の初期値は真値の近くにあるので最適化が安定します。「簡単だが効率は劣るMoM」と「効率的だが計算が要るMLE」を 組み合わせて 使う、実務的なベストプラクティスです。
最後に、同じデータに複数の分布をフィットして、形を比べてみましょう。
import matplotlib, matplotlib.pyplot as plt
for cand in ["Hiragino Sans", "Yu Gothic", "Noto Sans CJK JP", "IPAexGothic", "Meiryo"]:
if any(cand == f.name for f in matplotlib.font_manager.fontManager.ttflist):
plt.rcParams["font.family"] = cand
break
plt.rcParams["axes.unicode_minus"] = False
rng = np.random.default_rng(8)
x = rng.gamma(2.5, 1.6, 600) # 真はガンマ
m1, v = x.mean(), x.var()
# 3分布をMoMでフィット
g = mom_gamma(x)
gam = stats.gamma(g["k"], scale=g["theta"])
norm = stats.norm(m1, np.sqrt(v))
# 対数正規のMoM: s²=ln(1+v/m²), μ=ln(m)-s²/2
s2 = np.log(1 + v / m1**2)
ln = stats.lognorm(s=np.sqrt(s2), scale=np.exp(np.log(m1) - s2 / 2))
xs = np.linspace(0.001, x.max() * 1.05, 500)
plt.figure(figsize=(10, 5.4))
plt.hist(x, bins=40, density=True, color="#888", alpha=0.35,
edgecolor="white", label="標本（真はガンマ）")
plt.plot(xs, gam.pdf(xs), lw=2.4, color="#ff9900", label="ガンマ MoM")
plt.plot(xs, ln.pdf(xs), lw=2.0, ls="--", color="#1f77b4", label="対数正規 MoM")
plt.plot(xs, norm.pdf(xs), lw=2.0, ls=":", color="#2ca02c", label="正規 MoM")
plt.xlabel("x"); plt.ylabel("密度"); plt.legend()
plt.title("複数分布のMoMフィット比較 — 同じ平均・分散でも形が違う")
plt.grid(True, alpha=0.3); plt.tight_layout()
plt.show()
このグラフから重要な教訓が読み取れます。3つの分布はすべて 同じ平均・分散 をもつようにMoMでフィットしていますが、曲線の形はまったく違います。ガンマ（オレンジ）はヒストグラムにぴったり、対数正規（青破線）も似た非対称形でそこそこフィット、正規（緑点線）は対称なので負の領域にはみ出し、右の裾を捉えられていません。つまり「平均と分散を合わせる」だけでは分布の形まで一致する保証はなく、どの分布族を選ぶか が本質的に重要だということです。モーメント法は分布族を決めた上での当てはめ手法であり、分布族の選択そのものは別途データの形状（裾の重さ、台の範囲）を見て判断する必要があります。
まとめ
本記事では、モーメント法（積率推定）を直感から応用まで通して解説しました。
- 核心アイデア: 標本モーメント $\hat{m}_k$ を理論モーメント $\mu_k(\theta)$ に等置し、連立方程式を解く。パラメータ $p$ 個には $p$ 本のモーメントを使う
- 具体例: 正規は「平均と標本分散」、ガンマは $\hat{k} = \bar{x}^2/\hat{\sigma}^2$・$\hat{\theta} = \hat{\sigma}^2/\bar{x}$、ベータは $S = m(1-m)/v – 1$ から、一様は $m \pm \sqrt{3}\,s$ と、いずれも閉じた式で求まる
- MLEとの比較: モーメント法は簡単・高速だが効率はMLEに劣る（推定のばらつきが大きい）。一方でMLEの初期値として優秀
- 漸近性質: 大数の法則から一致性、中心極限定理＋デルタ法から漸近正規性をもつ
- GMM: モーメント条件がパラメータ数より多い過剰識別では、二次形式 $\bar{\bm{g}}^{\top}\bm{W}\bar{\bm{g}}$ を最小化する。尤度を仮定せず推定でき、計量経済学の標準ツールになっている
- 注意点: 一様分布のように台の境界がパラメータの分布では機能しにくい。また「平均・分散を合わせる」だけでは分布の形は一致しないので、分布族の選択が重要
モーメント法は、推定論への入り口として最適です。ここで学んだ「データの統計量を理論値に合わせる」という発想は、最尤推定やベイズ推定とも地続きであり、GMMを通じて計量経済学の最前線にもつながっています。
次のステップとして、以下の記事も参考にしてください。