【Rで統計入門】(3) モーメント法を使った確率分布のパラメタ推定

　確率分布のパラメータを推定する方法のひとつとして、カール・ピアソン（Karl Pearson, 1857–1936）が導入したモーメント法があります。直感的で計算がしやすいこの手法は、統計学の黎明期に広く用いられました。しかし、聡明な皆さんの中には、「最尤法に比べると単純すぎるのでは？」と感じる方もいるかもしれません。それでも、モーメント法は今でも現役で活躍する場面があり、特に、時系列のフラクタル解析や長時間相関解析などに応用されることもあります。私は、非ガウス性の解析にモーメント法を使う方法を提案しています。

　また、モーメント法を理解することは、より洗練された推定法である最尤推定法 (最尤法)を学ぶうえでも役立ちます。そして、何より興味深いのは、この推定法をめぐって繰り広げられたピアソンとロナルド・フィッシャー（Ronald A. Fisher, 1890–1962）の激しい対立です。二人の間には統計学の覇権をめぐる因縁があり、その物語は単なる数学的議論を超えて、まるでドラマのように展開していきます。今回は、モーメント法の基本を紹介するとともに、その歴史的背景として、ピアソンとフィッシャーの確執についてもお話ししたいと思います。数式を見たくない方は、「統計学をめぐる対立の物語：ピアソンとフィッシャー」まで、スキップしてください。それと、「母数」の意味についても触れておきました。

確率分布のパラメタ
モーメントとは
モーメント法による確率分布のパラメタ推定
統計学をめぐる対立の物語：ピアソンとフィッシャー
付録1. 「母数」という言葉の意味の変化について
付録2. ポアソン分布の１次モーメント

確率分布のパラメタ

　以前の記事で、統計の役割は、データを要約することと説明しましたが、観測データの平均や標準偏差を計算するだけでは、データについて何か分かったという気持ちにはなれません。なぜなら、そのデータは単なる数の集まりではなく、何らかの法則 (確率的な法則) に従って生じたかもしれないからです。そのような法則を発見できれば、同じような状況で将来何が起こるのかを予測したり、適切な対処を考えたりすることができます。そのようなデータを生み出す偶然の法則を記述するものとして確率分布があります。

　例えば、次のような現象について、それが従う確率分布のモデルが提案されています。

一定時間内に何回電話がかかってくるかについては、もし電話のタイミングがランダムで、１時間当たりに電話がかかってくる回数の平均が一定であれば、電話がかかって来る回数の分布は、ポアソン分布に従います。
機械や電子部品の寿命は、指数分布に従うことがあります。
試験の総合得点など、データが平均値の周りに分布する場合には、正規分布に従うことが多いです。

　もし、データの背後に、上で紹介したような確率分布があるのであれば、観測された標本から、確率分布の全体像をビシっと決定するのが、統計の役割です。

　ポアソン分布、指数分布、正規分布には、パラメタ (parameter: 母数)と呼ばれる未知数が含まれています。

ポアソン分布（Poisson distribution）

　ポアソン分布は、単位時間または単位空間内で発生する事象の回数を表す確率分布です。確率質量関数 (Probability Mass Function) $\displaystyle{ P(X = k) }$ は次のように表されます：

$\displaystyle{ P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0,1,2,\dots }$

ここで、パラメタ $\displaystyle{ \lambda \gt 0 }$ は、単位時間または単位空間あたりの平均発生回数（到着率）を表しています。

指数分布（Exponential distribution）

　指数分布は、独立した事象が発生するまでの時間を表す確率分布です。確率密度関数 (Probability Density Function: PDF) $\displaystyle{ f(x) }$ は次のように定義されます：

$\displaystyle{ f(x) = \lambda e^{-\lambda x}, \quad x \geq 0 }$

ここで、パラメタ (\lambda > 0)は、事象の発生率（単位時間あたりの平均発生率の逆数）を表しています。指数分布はポアソン過程の待ち時間分布としても知られています。

正規分布（Normal distribution）

　正規分布は、統計ではよく見るやつです。確率密度関数 (PDF)は次のように表されます：

$\displaystyle{ f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(x - \mu)^2}{2\sigma^2} \right), \quad -\infty < x < \infty }$

ここで、パラメタは、平均を表す $\displaystyle{ \mu }$ と、分散を表す $\displaystyle{ \sigma^ 2 }$ です。正規分布は、中心極限定理で説明されるように、多くの確率変数の合計や平均が従う分布として知られています。

モーメントとは

　モーメント (moment)は、確率変数の値を何乗かしたものの期待値 (平均)を求めたものです。観測されたデータを使った計算では、値を何乗かした場合の、平均を計算します。そして、 $q$ 乗した場合を、 $q$ 次モーメントと呼びます。 $q$ の値を変えることで、小さい値を強調したり、大きい値を強調したりして、分布の構造、あるいは、データの構造を調べることができます。

　具体的には、確率変数 $X$ の $q$ 次のモーメント $m _ q$ は次のように表されます。

$\displaystyle{ m_q = E[X^q] = \int_{-\infty}^{\infty} \! x^q \, f(x) \, dx }$

ここでは、 $f(x)$ は、確率変数 $X$ が従う確率密度関数です。

　観測された値

$\displaystyle{ \left\{x_1, x_2, \cdots, x_N \right\} }$

からモーメントを計算 (推定)したいときは、

$\displaystyle{ \hat{m}_q = \frac{1}{N} \sum_{i = 1}^N x_i^q }$

を計算します。観測されたデータから計算されたモーメントは、真の値ではなく推定なので、そのことを示すために、変数にハット"＾"を付けています。

　また、中心モーメント (central moment)というものもあり、これは、平均を中心としてモーメントを計算します。この場合、平均は、

$\displaystyle{ m_1 = E[X] }$

ですので、 $q$ 次モーメントの定義は、

$\displaystyle{ \mu_q = E \left[\left( X - m_1 \right)^q \right] }$

です。

　計算式からすぐにわかるように、2次中心モーメント $\displaystyle{ \mu _ 2 }$ は、分散と一致します。また、

$\displaystyle{ \mu_2 = m_2 - m_1^2 }$

が成り立ちます。

モーメント法による確率分布のパラメタ推定

　モーメント法を使って確率分布のパラメータを推定してみます。以下では、観測された $N$ 個のデータを

$\displaystyle{ \left\{x_1, x_2, \cdots, x_N \right\} }$

とします。

ポアソン分布のパラメタをモーメント法で求める

　ポアソン分布

$\displaystyle{ P(X=k)=\frac{\lambda^k e^{-\lambda}}{k!} }$

の1次モーメント $m _ 1$ (平均値)は、

$\displaystyle{ E[X]=\lambda }$

です (最後に付録でつけました)。

　したがって、標本平均を使ってパラメタ $\lambda$ を推定できます。つまり、

$\displaystyle{ \hat{\lambda}=\frac{1}{N} \sum_{i=1}^N x_i }$

です。観測されたデータから計算されたパラメタは、真の値ではなく推定なので、そのことを示すために、変数にハット"＾"を付けています。

指数分布のパラメータをモーメント法で求める

　指数分布

$\displaystyle{ f(x)=\lambda e^{-\lambda x} }$

の1次モーメント $m _ 1$ (平均値)は、

$\displaystyle{ E[X] = \int_0^\infty x \lambda e^{-\lambda x} dx = \frac{1}{\lambda} }$

です。

　したがって、標本平均を用いてパラメタを推定できます。

$\displaystyle{ \frac{1}{\hat{\lambda}} = \frac{1}{N} \sum_{i=1}^N x_i }$

正規分布のパラメタをモーメント法で求める

　ここまで書いてきて、当たり前の簡単な例しかあげられていないと感じましたが、それでもこのまま押し切ります。

　正規分布

$\displaystyle{ f(x) = \frac{1}{\sqrt{2\pi\sigma^2}} \exp \left( -\frac{(x - \mu)^2}{2\sigma^2} \right) }$

のパラメタの推定については、各パラメタの意味を知っていれば当たり前の関係式を使います。

　ということで、 $\mu$ は平均値 (1次モーメント)ですので、標本平均を使って計算できます。

　また、 $\sigma-2$ は分散 (2次中心モーメント)ですので、分散を計算すれば良いわけです。

　しかし、分散の計算って、標本分散なの？不偏分散なの？といろいろと計算式があってややこしいです。この点について、近いうちに説明します。

Rでパラメタ推定

以下は、各分布ごとにシミュレーションを行い、モーメント法によるパラメータ推定を実演するRスクリプトの例です。

# サンプルサイズを定義します
N <- 1000

# ---------------------------
# 1. ポアソン分布の場合
# ---------------------------
# 真のパラメタ λ を 5 とします
lambda_poisson <- 5

# ポアソン分布に従うデータを N 個生成します
data_poisson <- rpois(N, lambda = lambda_poisson)

# モーメント法では、1次モーメント（平均値）が λ に等しいので、標本平均を計算します
est_lambda_poisson <- mean(data_poisson)

cat("【ポアソン分布】\n")
cat("真のパラメタ λ =", lambda_poisson, "\n")
cat("モーメント法による推定値 λ-hat =", est_lambda_poisson, "\n\n")

# ---------------------------
# 2. 指数分布の場合
# ---------------------------
# 真のパラメタ λ を 2 とします
lambda_exp <- 2

# 指数分布に従うデータを N 個生成します (rate = λ)
data_exp <- rexp(N, rate = lambda_exp)

# 指数分布の1次モーメントは E[X] = 1/λ ですので、標本平均から λ を推定します
sample_mean_exp <- mean(data_exp)
est_lambda_exp <- 1 / sample_mean_exp

cat("【指数分布】\n")
cat("真のパラメタ λ =", lambda_exp, "\n")
cat("モーメント法による推定値 λ-hat =", est_lambda_exp, "\n\n")

# ---------------------------
# 3. 正規分布の場合
# ---------------------------
# 真のパラメタとして、平均 μ = 10、標準偏差 σ = 3 (分散 σ² = 9) とします
mu_normal <- 10
sigma_normal <- 3

# 正規分布に従うデータを N 個生成します
data_normal <- rnorm(N, mean = mu_normal, sd = sigma_normal)

# 正規分布では、1次モーメントが平均 μ、2次中心モーメントが分散 σ² となるため、
# 標本平均と標本分散を計算して推定します
est_mu_normal <- mean(data_normal)
est_sigma2_normal <- var(data_normal)  # Rのvar関数は不偏分散 (n-1で割る) を返しますが、ここでは推定の例として使用します

cat("【正規分布】\n")
cat("真のパラメタ μ =", mu_normal, ", σ² =", sigma_normal^2, "\n")
cat("モーメント法による推定値 μ-hat =", est_mu_normal, ", σ²-hat =", est_sigma2_normal, "\n")

スクリプトの内容についての説明

　上のRスクリプトは、モーメント法を用いて各確率分布のパラメータを推定する数値実験を行っています。まず、サンプルサイズ(N=1000)を設定し、結果の再現性を確保しています。

ポアソン分布の場合
真のパラメタ( \lambda )を5に設定し、rpois関数でポアソン分布に従うデータを生成します。モーメント法では、ポアソン分布の平均が( \lambda )と等しいため、生成したデータの標本平均を計算し、推定値として出力します。
指数分布の場合
真のパラメタ( \lambda )を2に設定し、rexp関数で指数分布のデータを生成します。指数分布では1次モーメント（平均値）は(1/\lambda)であるため、標本平均の逆数をとることで( \lambda )の推定値を求めます。
正規分布の場合
真の平均( \mu )を10、標準偏差( \sigma )を3（分散( \sigma²=9 )）に設定し、rnorm関数で正規分布のデータを生成します。正規分布では、平均は1次モーメント、分散は2次中心モーメントに対応するため、標本平均と標本分散をそれぞれ計算して、( \mu )と( \sigma² )を推定します。

この実験を通じて、モーメント法がどのようにしてデータからパラメータを推定するか、それを体験してください。

統計学をめぐる対立の物語：ピアソンとフィッシャー

　カール・ピアソンは、自らを統計学の帝王と信じて疑いませんでした。19世紀末から20世紀初頭にかけて、彼は統計学を「科学の言語」として確立し、多くの弟子たちを従え、イギリスの統計学界を支配していました。彼の提唱したモーメント法は、データの平均や分散などの基本的な統計量を利用して確率分布のパラメータを推定する便利な手法であり、当時の統計学者たちの間で広く使われていました。

　「データの特徴をそのまま使ってパラメータを推定すればいい。計算もシンプルだし、十分に役に立つだろう」とピアソンは考えていました。彼にとって、この方法こそが統計学を実用的なものにする道筋だったのです。モーメント法は多くの場面でそれなりにうまく機能し、ピアソンは満足していました。しかし、そんな彼の王国に、一人の若き革命家が現れることになりました。

　ロナルド・フィッシャーは、ピアソンより30年以上若い世代の数学者でした。彼は鋭い頭脳を持ち、数学的厳密さにこだわる人物でした。ある日、彼はピアソンのモーメント法を眺めながら、ふと思いました。

　「これは……あまりに雑すぎる」

　フィッシャーの目には、モーメント法が単なる近似的な方法にしか見えませんでした。確かに計算は簡単ですが、果たしてそれが本当に最も良い推定なのか、という疑問が湧いてきたのです。
「推定量というのは、もっと精密で、理論的に最もよい方法であるべきだ」とフィッシャーは考えました。そして彼は、まったく新しい手法を生み出しました。それが「最尤法（Maximum Likelihood Estimation）」でした。

　フィッシャーのアイデアは、「データが観測される確率を最大にするようなパラメータを求めるべきだ」というものでした。直感的に聞くともっともらしいですが、当時の統計学の世界ではこのようなアプローチは斬新でした。最尤法は、データが生じる確率（尤度）を最大にするパラメータを数学的に求める手法であり、理論的に美しく、計算の精度も高かったのです。

　「これこそが真の統計学だ！」とフィッシャーは意気揚々としていました。彼はこの手法を論文にまとめ、世間に発表しました。しかし、ここで彼の前に立ちはだかったのが、ピアソンでした。

　ピアソンは、フィッシャーの論文を目にすると、鼻で笑いました。「最尤法だと？そんな複雑なことをしなくても、モーメント法で十分じゃないか。実用性が大事なんだよ」。彼は、フィッシャーの論文が掲載されるのを妨害し、学界での発表の場を与えようとしませんでした。

　フィッシャーはプンプンに憤りました。「じじいはただ、自分の地位を守りたいだけじゃないか」

　フィッシャーはそれでも諦めませんでした。彼は数学的に最尤法が「一致性」「効率性」「不偏性」という良い性質を持つことを示し、最尤法こそが最も信頼できる推定法であることを証明していきました。

　ピアソンも負けてはいませんでした。彼は自分の統計学的王国を守るため、フィッシャーを「無礼な若造」とみなし、彼のアイデアをことごとく否定し続けました。論文の掲載を妨害したり、彼の研究成果を公然と批判したりしたのです。ピアソンの弟子たちも師の意向を汲み取り、フィッシャーの最尤法を「理論的には面白いが、実際には使えない」と切り捨てました。

　しかし、時間が経つにつれて、学界はフィッシャーの考えに傾き始めました。最尤法の理論的な美しさと精度の高さは、多くの数学者や統計学者を魅了し、少しずつ支持を集めていったのです。そして最終的に、最尤法は統計学の主流となり、モーメント法は「計算が楽な代替手段」として二番手に追いやられることになりました。

　ピアソンとフィッシャーの確執は生涯続きました。ピアソンが亡くなるまで、二人は互いを認めることはありませんでした。しかし、皮肉なことに、今日の統計学ではピアソンもフィッシャーもどちらも偉大な業績を残したと評価されているのです。ピアソンは統計学を体系化し、多くの基礎的な概念を築きました。そしてフィッシャーは、統計学をより精密で理論的なものへと進化させました。

　歴史の勝者はフィッシャーだったかもしれません。しかし、もしピアソンがいなければ、統計学はここまで発展しなかったでしょう。
「統計学の父」と呼ばれる二人の男。彼らの確執こそが、現代の統計学を形作ったのです。

付録1. 「母数」という言葉の意味の変化について

　統計学において「母数（parameter）」は、もともと母集団全体の特性を表す数値のことを指します。例えば、ある国の全人口の平均身長や、すべての工場で作られた製品の平均寿命といった値が「母数」に当たります。これは、私たちが調査を通じて知りたい真の値であり、現実には直接測定するのが難しいものです。そのため、通常は一部のデータ（サンプル）をもとに推定することになります。

　しかし、最近では「母数」という言葉が「サンプルサイズ（sample size）」、つまり調査対象の数という意味で使われることが増えてきています。たとえば、あるアンケート調査で「今回の母数は500人です」といった表現を見かけることがあります。この使い方は、本来の統計学的な意味とは異なり、「調査対象の人数」という意味合いで使われています。

　なぜこのような意味の変化が起こったのでしょうか？一つの理由として、「母数」という言葉が「データの規模」や「対象の数」といった漠然とした意味で受け取られやすいことが挙げられます。統計を専門的に学んでいない人にとって、「母数」という言葉は「母集団の特性」という抽象的な概念よりも、「対象の数」のような具体的なイメージのほうがしっくりくるのかもしれません。また、実際の調査現場では、サンプルサイズの規模が重要なポイントとなるため、その影響も考えられます。

　このような意味の変化に対して、統計学の専門家の中には違和感を持つ人も少なくありません。学術的には「母数＝母集団の特性」「サンプルサイズ＝調査対象の数」と明確に区別されており、誤用が広まることで誤解が生じる可能性もあります。一方で、言葉は時代とともに変化するものでもあります。現場で使われるうちに、新しい意味が定着していくことも珍しくありません。

「母数」という言葉の変化は、統計学の専門的な概念と一般社会での実用的な使われ方の違いを象徴する一例と言えるでしょう。正確な表現を心がけることは大切ですが、一方で、言葉の変遷が生まれる背景を理解することもまた、興味深い視点かもしれません。

付録2. ポアソン分布の１次モーメント

　ポアソン分布の１次モーメント (期待値)の計算過程を示しておきます。

$\displaystyle{ \begin{align} E[X] & =\sum_{k=0}^{\infty} k\, P(X=k)\\ &= \sum_{k=0}^{\infty} k \, \frac{\lambda^k e^{-\lambda}}{{k}!} \\ &= 0 + \sum_{k=1}^{\infty} k \, \frac{\lambda^k e^{-\lambda}}{{k}!} \qquad(0! = 1) \\ &= \sum_{k=1}^{\infty} \frac{\lambda^k e^{-\lambda}}{(k-1)!} \\ &= \sum_{k=1}^{\infty} \frac{\lambda \cdot \lambda^{k-1} e^{-\lambda}}{(k-1)!} \\ &= \lambda \sum_{k=1}^{\infty} \frac{\lambda^{k-1} e^{-\lambda}}{(k-1)!} \\ &= \lambda \sum_{j=0}^{\infty} \frac{\lambda^{j} e^{-\lambda}}{j!} \qquad (k-1 = j)\\ & = \lambda \cdot 1 \\ & = \lambda \end{align} }$

　確率はすべて足すと1ですので、

$\displaystyle{ \begin{align} \sum_{k=0}^{\infty} P(X=k) = \sum_{k=0}^{\infty} \frac{\lambda^{k} e^{-\lambda}}{k!} = 1 \end{align} }$

を最後の変形で使いました。また、この式は、以下のように計算すれば求められます。

$\displaystyle{ \begin{align} \sum_{k=0}^{\infty} P(X=k) &= \sum_{k=0}^{\infty} \frac{\lambda^{k} e^{-\lambda}}{k!}\\ &= e^{-\lambda} \sum_{k=0}^{\infty} \frac{\lambda^{k} }{k!} \\ &= e^{-\lambda} \left\{1 + \frac{\lambda^1}{1!} + \frac{\lambda^2}{2!} + \cdots \right\} \\ &= e^{-\lambda} \cdot e^{\lambda}\\ &= 1 \end{align} }$

ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります