2019-03-15

クーポンコレクター問題の確率分布を解き明かす

確率論

クーポンコレクター問題というものをご存知だろうか？これは、例えば6種類のおもちゃが出るガシャポン *1があったとして、何回くらい引けば全種類引き当てる事ができるか？というようなことを考える問題である。

この問題に対して、平均や分散がどうなるかということは非常によく語られることである[1]。また、時として不等式による評価について議論している記事を見かけることもある[2]。

しかし、肝心の確率分布については議論される事が非常に少ない。あまりにも記事を見かけないので、私は当初クーポンコレクター問題の確率分布を求めることは不可能なのではないかとさえ思っていた。

それでもめげずに調べ続けた結果、私はついに確率分布について結論を出しているページを見つけた。本稿ではそれを紹介し、喜びを分かち合いたいと思う。

第2種スターリング数

クーポンコレクター問題の確率分布を求めるためには、第2種スターリング数について理解しておく必要がある。少し長いが、定義をWikipediaから引用する[3]。

第2種スターリング数

第2種スターリング数 (Stirling number of the second kind) $\{{\textstyle {n \atop k}}\}$ は、 $x^{n}$ を下降階乗冪 $x^{\underline {k}}\equiv x\,(x-1)(x-2)\cdots (x-k+1)$ の級数:

$\displaystyle{ x^{n}=\sum _{k=0}^{n}\left\{{n \atop k}\right\}\,x^{\underline {k}} }$

で展開したときの展開係数として定義される。この定義では、 $0\leq k\leq n$ である。便宜上、 $\{{\textstyle {0 \atop 0}}\}=1$ と定義する。第2種スターリング数は

$\displaystyle{ \left\{{n \atop k}\right\}=\left\{{n-1 \atop k-1}\right\}+k\,\left\{{n-1 \atop k}\right\} }$

なる漸化式で計算できる。

定義も大切なのだが、第2種スターリング数は以下のように特徴付けられるということが重要である[3]。

第2種スターリング数の特徴付け

第2種スターリング数 $\{{\textstyle{n \atop k}}\}$ は、組合せ数学において、番号づけされた $n$ 個の要素をグループ $k$ 個に分割する組み合わせの数を与える。分割する要素は番号付けされているので個別に区別できるが、グループは順序を特に区別しないものとする。

クーポンコレクター問題の確率分布

準備が整ったので、確率分布を求めてみよう。以下の議論は全て[4][5]を参考にした。

集めるクーポンの種類を $k$ 種類とし、初めて全種類のクーポンを取得できるまでの試行回数を表す確率変数を $N$ とする。この時、 $n$ 回目の試行で初めて全種類のクーポンを取得できる確率 $P(N = n)$ を求める。

$n$ 回の試行により出現し得るクーポンの出方の総数は $k^n$ である。もし $n$ 回目の試行で初めて全種類のクーポンが出るパターンの総数が分かれば、それを $k^n$ で割ったものが求める確率である。

$n$ 回目の試行で初めて全種類のクーポンが出揃うということは、 $n-1$ 回目の時点で $k-1$ 種類のクーポンがすでに出ている必要がある。このパターンの総数を求めるために、 $n-1$ 回目までの試行に対して1から順に番号を割り振る。そして、これらの番号付けられた試行を $k-1$ 個のグループに分ける。同じグループに分けられた試行については同じクーポンが得られたと考える。この分け方のパターンの総数は第2種スターリング数 $\{\textstyle{n-1 \atop k-1}\}$ となる。

第2種スターリング数は分けられたグループの間の順序は区別しないが、今はクーポンの種類は区別されるので、 $(k-1)!$ を掛ける。さらに、 $n-1$ 回目までに出る $k-1$ 種類のクーポンは $k$ 種類のうちどれであるかは問わないので、全てのパターンをカウントする必要がある。そのため、 $\left(\textstyle{k \atop k-1}\right) = k$ を掛ける。結局、 $n-1$ 回目の時点で $k-1$ 種類のクーポンが出現するパターンの総数は $k! \{\textstyle{n-1 \atop k-1}\}$ となる。

$n$ 回目の試行でまだ出ていない最後のクーポンが出るパターンの総数は1なので、 $n$ 回目の試行で初めて全種類のクーポンが出るパターンの総数は $k! \{\textstyle{n-1 \atop k-1}\}$ となる。以上の議論により次の式を得る。

$\displaystyle{ P(N = n) = \frac{k!}{k^n} \{\textstyle{n-1 \atop k-1}\} }$

計算してみよう

$n, k$ にいくつか具体的な数値を入れて $P(N = n)$ を計算してみよう。手計算ではやっていられないので、pythonを利用する。第2種スターリング数はsympyのstirling関数を用いて計算することが出来るので、それを利用する[6]。数値だけ示しても分かり辛いので、結果をmatplotlibでグラフ化した。Python自体はpydroid3を使って実行した[7]。なお、sympyなどの必要なパッケージはあらかじめpipでインストールしておく必要がある。

計算に用いたソースコードを以下に示す。

from sympy.functions.combinatorial.numbers import stirling
import matplotlib.pyplot as plt
import math

startCouponKind = 2
numCouponKind = 6
numTrial = 30
for k in range(startCouponKind, numCouponKind+1):
  x = []
  y = []
  for n in range(1, numTrial+1):
    stir = stirling(n-1, k-1, kind = 2)
    prob = math.factorial(k) * stir / (k**n)
    x.append(n)
    y.append(prob)
  plt.plot(x, y, label='k = {}'.format(k))
plt.legend()
plt.savefig('/storage/emulated/0/coupon.png')

得られたグラフを以下に示す。

f:id:peng225:20190314084353j:plain

途中でピークを持ち、右側に裾野が広がっている様子が見て取れる。 $n$ が小さいときに確率が0となる領域があるが、これはクーポンの種類よりも試行回数が少ないときの様子が表れているものである。

また、 $k$ の増え方に対して、ピークの位置が右に移動していく速度の方がやや速いように見える。実際、全種類コンプリートするまでの試行回数の期待値は $O(k \mathrm{log} k)$ なので、これは理論的にも辻褄が合っている[1]。

まとめ

本稿ではクーポンコレクター問題の確率分布を明らかにし、実際に計算を行った結果を示した。クーポンコレクター問題ではパターンを数え上げる事が難しかったが、その難しさを既知の概念である第2種スターリング数に押し込めることで、理論的にスッキリとした結論を得ることができた。

これまで組み合わせ数学には興味がなかったが、本稿を書き上げるうちにその面白さを垣間見る事ができた。またいずれ体系的に勉強しよう。

参考

[1] クーポンコレクター問題 - Wikipedia
[2] コンプガチャの数理 -コンプに必要な期待回数の計算方法について- - doryokujin's blog
[3] スターリング数 - Wikipedia
[4] combinatorics - Probability distribution in the coupon collector's problem - Mathematics Stack Exchange
[5] CDF of probability distribution with replacement - Mathematics Stack Exchange
[6] Combinatorial — SymPy 1.3 documentation
[7] https://play.google.com/store/apps/details%3Fid%3Dru.iiec.pydroid3%26hl%3Dja%26referrer%3Dutm_source%253Dgoogle%2526utm_medium%253Dorganic%2526utm_term%253Dpydroid3%26pcampaignid%3DAPPU_1_5SOHXJe-E4-Pr7wPpLSfwAE

*1:これの呼び方は地方によって差があるような気がするが、自分の流儀で呼ばせて頂く。

2019-03-01

いくつかのLie群がLie群であることを定義に戻って確かめる

Lie群

Lie 群は難しい。この理由の1つは、議論の前提となる領域が広いことにあると思われる。Lie群とは群であり多様体であるような数学的対象である。そのため、定義を理解するだけで群論と多様体の知識が求められる。また、Lie群の教科書で最初に扱われるような基本的なLie群は行列群である。しかも、そのコンパクト性に着目した議論も多い。そのため、線形代数と位相空間の基礎的な事項も理解しておくことが望ましい。

繰り返すが、Lie群は難しい。私はここ最近Lie群を勉強し始めて、この事実を痛感している。こういう時は足元を一歩ずつ踏み固めて行くしかない。その一環として、本稿ではいくつかの基本的なLie群について、それらが本当にLie群になっていることを定義に照らし合わせて確認してみる。

本稿では私の独断で以下の2つのLie群を扱う。

一般線形群 $\mathrm{GL}(n, \mathbb{R})$
直交群 $\mathrm{O}(n)$

準備

Lie群の定義

Lie群の定義を[1]より引用する。

Lie群

多様体 $G$ が群構造を持ち、その群演算
$\displaystyle{ G \times G \to G;\ (x, y) \mapsto xy, \ \ G \to G;\ x \mapsto x^{-1} }$
が $C^{\infty}$ 級写像になるとき、 $G$ をLie群と呼ぶ。

多様体上の写像が $C^{\infty}$ 級であるということ

Lie群の定義の中で $C^{\infty}$ 級写像という言葉が出てきた。この定義を[2]より引用する。

1点においてC^s級

連続写像 $f: M \to N$ が, 1点 $p \in M$ において $C^s$ 級であるとは, $p$ を含む $M$ の $C^r$ 級座標近傍 $(U; x_1, \cdots , x_m)$ と $f(p)$ を含む $N$ の $C^r$ 級座標近傍 $(V; y_1, \cdots , y_n)$ が存在して,
(1) $f(U) \subset V$ かつ
(2) $(U; x_1, \cdots , x_m)$ と $(V; y_1, \cdots , y_n)$ に関する $f$ の局所座標表示が $C^s$ 級である,
この2つの条件がなりたつことである. ただし, $1 \le s \le r$ .

C^s級写像

$f: M \to N$ が $C^s$ 級写像 (mapping of class $C^s$ ) であるとは, $M$ の各点 $p$ において $f$ が $C^s$ 級であることである.

上記定義において $s = \infty$ とすれば $C^{\infty}$ 級写像の定義となる。

部分多様体に関する定理

後で使う定理について説明するために、正則点・臨界点、及び正則値・臨界値の定義を述べておく[2]。

正則点・臨界点

$M$ の点 $p$ における $f$ の微分
$\displaystyle{ (df)_p : T_p (M) \to T_{f(p)} (N) }$
が‘上へ’の線形写像であるとき, $p$ を $f$ の正則点 (regular point) とよぶ. $(df)_p$ が‘上へ’の線形写像でないとき, $p$ を $f$ の臨界点 (critical point) という.

正則値・臨界値

$f : M \to N$ の臨界点全部の集合 (臨界点集合) を $C_f$ で表す ( $C_f$ は $M$ の部分集合である) .
$C_f$ の像 $f(C_f)$ に属する $N$ の点 $q$ を $f$ の臨界値 (critical value) とよぶ. 臨界値でない $N$ の点を $f$ の正則値 (regular value) という.

以下の定理を後の議論で使用する[2]。

部分多様体に関する定理

点 $q \in N$ が $C^r$ 級写像 $f: M \to N$ の正則値で, $f^{-1}(q) \neq \phi$ であるとすると, 逆像 $f^{-1}(q)$ は $M$ の $(m-n)$ 次元 $C^r$ 級部分多様体である.

Lie群であることの確認

以上で準備が整ったので、Lie群であることの確認に移る。本稿では群であることの確認はサボり、それぞれのLie群について以下の3点を確認した。

$C^{\infty}$ 級多様体であること
群演算が $C^{\infty}$ 級写像であること
逆元を取る操作が $C^{\infty}$ 級写像であること

一般線形群 $\mathrm{GL}(n, \mathbb{R})$

$C^{\infty}$ 級多様体であること

$\mathbb{R}$ 上 $n$ 次正方行列全体の集合を $\mathrm{M}(n, \mathbb{R})$ と書く。 $\mathrm{GL}(n, \mathbb{R})$ は $\mathrm{M}(n, \mathbb{R})$ の部分集合のうち、以下のように表されるものである。

$\displaystyle{ \{A \in \mathrm{M}(n, \mathbb{R}) | \mathrm{det} (A) \neq 0\} }$

実は、 $\mathrm{GL}(n, \mathbb{R})$ は $\mathrm{M}(n, \mathbb{R})$ の開部分集合となる。詳細は[3]の命題1.17に譲り、ここでは概要だけ説明する。まず、 $\mathrm{det} : \mathrm{M}(n, \mathbb{R}) \to \mathbb{R}$ は連続写像である。このとき、 $\mathrm{GL}(n, \mathbb{R}) = \mathrm{det}^{-1}\left(\{x \in \mathbb{R} | x \ne 0\}\right)$ と書ける。 $\{x \in \mathbb{R} | x \ne 0\}$ は $\mathbb{R}$ の開集合なので、連続写像 $\mathrm{det}$ による逆像も開集合となる。

ここで、 $\mathrm{M}(n, \mathbb{R})$ に属する行列の各成分を座標と見なすと、これは $\mathbb{R}^{n^2}$ と同一視できる。そのため、
$\mathrm{M}(n, \mathbb{R})$ は $C^{\infty}$ 級多様体となる。その開部分集合も $C^{\infty}$ 級多様体となるので、 $\mathrm{GL}(n, \mathbb{R})$ は $C^{\infty}$ 級多様体となる。このような多様体を開部分多様体と呼ぶ[2]。

群演算が $C^{\infty}$ 級写像であること

$\mathrm{GL}(n, \mathbb{R})$ の群演算は通常の行列の掛け算である。この写像を $\mathrm{mul} : \mathrm{GL}(n, \mathbb{R}) \times \mathrm{GL}(n, \mathbb{R}) \to \mathrm{GL}(n, \mathbb{R})$ と表す。以下で $\mathrm{mul}$ が $C^{\infty}$ 級であることを示す。

$\mathrm{GL}(n, \mathbb{R})$ はそれ自身 $\mathbb{R}^{n^2}$ の開集合と見なせるため、そのアトラスとして自分自身だけを座標近傍として含む、要素数1の集合族を取れる。 $\mathrm{GL}(n, \mathbb{R}) \times \mathrm{GL}(n, \mathbb{R})$ についても、積多様体のアトラスの定義を考えると、やはりアトラスとして自分自身だけを要素として含む集合族を取れる。よって $\mathrm{GL}(n, \mathbb{R})$ 、 $\mathrm{GL}(n, \mathbb{R}) \times \mathrm{GL}(n, \mathbb{R})$ は共に $C^{\infty}$ 級多様体である。

任意の点 $(A, B) \in \mathrm{GL}(n, \mathbb{R}) \times \mathrm{GL}(n, \mathbb{R})$ について、 $\mathrm{mul}( (A, B) ) = AB$ の各成分は $A, B$ の成分の多項式となる。よって $\mathrm{mul}$ は $(A, B)$ において $C^{\infty}$ 級である。 $(A, B)$ は任意なので、 $\mathrm{mul}$ は $C^{\infty}$ 級写像である。

逆元を取る操作が $C^{\infty}$ 級写像であること

$\mathrm{GL}(n, \mathbb{R})$ の元の逆元を取る操作とは、逆行列を求める操作を意味する。この写像を $\mathrm{inv} : \mathrm{GL}(n, \mathbb{R}) \to \mathrm{GL}(n, \mathbb{R})$ と表す。以下で $\mathrm{inv}$ が $C^{\infty}$ 級であることを示す。

任意の点 $A \in \mathrm{GL}(n, \mathbb{R})$ について、 $\mathrm{inv}(A) = A^{-1}$ である。逆行列は余因子行列を行列式で割ることによって得られるため、 $A^{-1}$ の各成分は $A$ の成分の有理式となる。よって $\mathrm{inv}$ は $A$ において $C^{\infty}$ 級である。 $A$ は任意なので、 $\mathrm{inv}$ は $C^{\infty}$ 級写像である。

$\mathrm{inv}$ の定義域は $\mathrm{GL}(n, \mathbb{R})$ なので、行列式は常に0ではないことに注意されたい。

直交群 $\mathrm{O}(n)$

$C^{\infty}$ 級多様体であること

これを示すのは思いのほか難しいため、証明のアウトラインだけ述べることにする。詳細は[4]を参照されたい。

$n$ 次実対称行列全体の集合を $\mathrm{S}(n)$ とする。 $\mathrm{S}(n)$ は $\mathbb{R}^{n(n+1)/2}$ と同一視できるため、 $n(n+1)/2$ 次 $C^{\infty}$ 級多様体である。

この時、以下のような写像を考える。

$\displaystyle{ f : \mathrm{M}(n, \mathbb{R}) \to \mathrm{S}(n), A \mapsto AA^T }$

$f$ は $C^{\infty}$ 級写像である。 $A$ が直交行列のとき、その逆行列は $A^T$ となる。そのため、 $\mathrm{O}(n) = f^{-1}(I)$ となる。ここで、 $I$ は単位行列である。

もし $I$ が $f$ の正則値であれば、先ほど提示した定理により $\mathrm{O}(n)$ は $n(n-1)/2$ 次 $C^{\infty}$ 級多様体となる。そのためには、 $\mathrm{O}(n)$ の全ての点における微分が全射になれば良い。具体的に微分計算を行い、それが全射であることを確かめるアプローチになるが、体力の限界なのでこれ以降は[4]にお任せする。

群演算が $C^{\infty}$ 級写像であること

$\mathrm{O}(n)$ の群演算は $\mathrm{GL}(n, \mathbb{R})$ における群演算を $\mathrm{O}(n) \times \mathrm{O}(n)$ に制限した写像 $\mathrm{mul} |_{\mathrm{O}(n) \times \mathrm{O}(n)}$ によって与えられる。これは包含写像 $i : \mathrm{O}(n) \times \mathrm{O}(n) \to \mathrm{GL}(n, \mathbb{R}) \times \mathrm{GL}(n, \mathbb{R})$ と $\mathrm{mul}$ によって以下のように書ける。

$\displaystyle{ \mathrm{mul} |_{\mathrm{O}(n) \times \mathrm{O}(n)} = \mathrm{mul} \circ i }$

$i$ と $\mathrm{mul}$ はどちらも $C^{\infty}$ 級であるため、それらの合成写像も $C^{\infty}$ 級となる。よって $\mathrm{O}(n)$ の群演算は $C^{\infty}$ 級写像である。

逆元を取る操作が $C^{\infty}$ 級写像であること

群演算とほぼ同じように制限写像を考えれば良い。詳細は割愛する。

まとめ

本稿では一般線形群と直交群が共にLie群であることを確認した。どちらも非常に有名なLie群でありながら、それを示すのは意外と大変だった。Lie群を理解するまでの道のりは険しい。

参考

[1] http://www.math.tsukuba.ac.jp/~tasaki/lecture/ln2010/2010t.pdf
[2]

多様体の基礎 (基礎数学5)

作者: 松本幸夫
出版社/メーカー: 東京大学出版会
発売日: 1988/09/22
メディア: 単行本
購入: 7人クリック: 36回
この商品を含むブログ (33件) を見る

[3] http://www.geocities.jp/aoirei2002/math/papers/topogp.pdf
[4] http://math.uchicago.edu/~may/REU2014/REUPapers/Rouse.pdf

2019-01-27

今こそJordan標準形と向き合う

線形代数

線形代数を勉強して、Jordan標準形という言葉を耳にしたことがない人はいないだろう。Jordan標準形とは、ざっくり言えば行列の対角化を一般化したようなものである。行列の対角化はいつでもできるとは限らず、報われない行列たちが存在する。一方、Jordan標準形は常に存在することが知られており、まさに行列界の救世主と言える。

私もかつて大学院入試の際にJordan標準形について勉強したことはある。しかし、専攻が情報系というのもあって、当時勉強したのは主にJordan標準形への変形方法だけで、その背後にある数学的な面白さについては理解していなかった。

そこで、本稿ではJordan標準形とは何なのか、その理論的な詳細について考えてみる。

行列は線形写像の映し鏡

始めに、本稿を読み進めるにあたって重要となる考え方について述べておく。

ベクトル空間上の線形写像は行列によって表現することができる。特に、線形変換は線形写像であるから、これも行列によって表現可能である。対角化やJordan標準形を考える上では、この逆を考えることが重要である。すなわち、ある行列に対して、それが何らかの線形変換の表現行列なのだと捉えるのである。

行列というのは線形変換をある基底に対して表現したものに過ぎず、基底の取り方によって姿を変え得る。唯一不変なのはその背後にある線形変換であり、線形変換を通して行列の振る舞いを考えることがJordan標準形の理解へと繋がる。

線形変換を通したベクトル空間の分解

固有空間への分解

$n$ 次元ベクトル空間 $V$ に対して、ある線形変換 $T$ が定められており、 $m$ 個の固有値 $\lambda_1, \lambda_2, \cdots , \lambda_m$ を持つとする。この時、各固有値について、それに対応する固有空間 $W_1, W_2, \cdots, W_n$ が存在する。ここでは線形変換の固有値、及び固有空間について述べているのであって、行列については一切触れていないことに注意されたい。

固有空間の基底は固有ベクトルであるから、固有空間に属する任意のベクトルに $T$ を作用させると、それらは単に元のベクトルの固有値倍となる。ゆえに、固有空間は $T$ 不変な部分空間となる。

全ての固有空間の次元の和が $n$ に等しい場合、以下の式が成立する。

$\displaystyle{ \begin{eqnarray} && V = W_1 \oplus W_2 \oplus \cdots \oplus W_n \\ && T(W_i) \subseteq W_i \end{eqnarray} }$

広義固有空間への分解

さて、いつもこのようになっていれば良いのだが、残念ながらそうはならないケースが存在する。すなわち、全ての固有空間の次元の和が $n$ を下回るような場合である。このような場合でも、何とかうまく $V$ を直和分解出来ないだろうか？

実は、これは可能である。 $W_i$ が固有空間であるという条件を緩和して、 $T$ 不変性と直和分解されるという性質だけを担保することで、 $V$ を先ほどと全く同じような直和分解の形に持ち込むことができる。

では、具体的にどのような部分空間に分解してやれば良いだろうか？それを考えるために、まずは固有ベクトル ${\bf x} \in W_i$ について、以下のような式変形を行う。

$\displaystyle{ \begin{eqnarray} T{\bf x} &=& \lambda_i {\bf x} \\ (T - \lambda_i I){\bf x} &=& {\bf 0} \end{eqnarray} }$

これより、固有ベクトルというのは $T - \lambda_i I$ を作用させると零ベクトルになるようなベクトルであると言える。または、簡潔に $W_i = \mathrm{Ker} (T - \lambda_i I)$ である。

ここから着想を得て、 $T - \lambda_i I$ だけでなく、 $(T - \lambda_i I)^k$ ( $k$ は任意の自然数) を作用させて零ベクトルになるようなベクトルの集合を考えてみる。これが実は $V$ の部分空間となっており、さらに $T$ 不変性と直和分解されるという性質を満たしている。このようにして得られる空間を広義固有空間と呼ぶ。これを $W_i'$ と表記すると、やはり簡潔に $W_i' = \mathrm{Ker} \{(T - \lambda_i I)^k\}$ と書ける。

少し確認してみよう。まず $T$ 不変性について、以下の式を考える。

$\displaystyle{ \begin{eqnarray} (T - \lambda_i I)^k (T{\bf x}) &=& (T - \lambda_i I)^{k-1}(T^2 - \lambda_i T) {\bf x} \\ &=& (T - \lambda_i I)^{k-1}T(T - \lambda_i) {\bf x} \\ &=& \cdots \\ &=& T(T - \lambda_i I)^k {\bf x} \\ &=& {\bf 0} \end{eqnarray} }$

$(T - \lambda_i I)^k$ を作用させて零ベクトルになったので、 $T{\bf x} \in W_i'$ が言える。

続いて $V$ が各 $W_i'$ の直和になることについてだが、これは思ったより証明が大変なので、ここではサボって[1]に譲ることにする。

直和分解について一点だけ注意事項を述べておく。ここまで広義固有空間への直和分解について説明したが、実は広義固有空間はさらに直和分解できる場合がある。具体的には、後述するJordan鎖の数だけさらに直和分解可能である。詳細はここでは述べないが、[2]などが参考になるだろう。

広義固有空間の基底

次に、 $W_i'$ の基底について考えてみよう。 $W_i'$ の次元を $m_i$ とする。固有ベクトル ${\bf x}_1 \in W_i'$ について、 $(T - \lambda_i I) {\bf x}_2 = {\bf x}_1$ を満たすベクトル ${\bf x}_2$ を考える。このような ${\bf x}_2$ は存在するかもしれないし、しないかもしれない。もし存在すれば、 $(T - \lambda_i I)^2 {\bf x}_2 = {\bf 0}$ となるため、 ${\bf x}_2 \in W_i'$ が言える。しかも、これは ${\bf x}_1$ のスカラー倍でもない。そのため、 $\{{\bf x}_1, {\bf x}_2\}$ は一次独立である。

帰納的に、 $(T - \lambda_i I) {\bf x}_{j+1} = {\bf x}_j$ となるベクトル ${\bf x}_{j+1}$ を考える。すると、 $(T - \lambda_i I)^{j+1} {\bf x}_{j+1} = {\bf 0}$ となるため ${\bf x}_{j+1} \in W_i'$ が言える。しかも、 $\{{\bf x}_1, {\bf x}_2\, \cdots , {\bf x}_{j+1}\}$ は一次独立になる。

基底は最大で $m_i$ 個しか取れないので、この操作はどこかで頭打ちとなる。このようにして得られるベクトルの列をJordan鎖という。Jordan鎖は1つ以上存在し、全てのJordan鎖を構成する全ベクトルを寄せ集めると、これは $W_i'$ の基底となる。これを広義固有ベクトルと呼ぶ。

広義固有空間の構造を完全に明らかにするためには、Jordan鎖は長さいくつのものが何本存在するのかを知る必要がある。これは $\mathrm{Ker} \{(T - \lambda_i I)\}, \mathrm{Ker} \{(T - \lambda_i I)^2\}, \cdots$ の次元を順に計算し、それらがどのように増えていくかを調べれば分かる。が、込み入った話になるので詳細は[3]を参照されたい。

広義固有ベクトルを基底とした場合の表現行列

以上、線形変換の性質についていろいろと述べたが、ここからいよいよJordan標準形の話に入っていく。

いきなり天下り的だが、線形変換 $T$ について、全ての固有値に対する広義固有ベクトルを全て集め、それらを基底とした場合の表現行列について考えてみよう。同じ広義固有空間から抽出した広義固有ベクトルは隣り合うように並べて、その中でさらに同じJordan鎖に属するベクトルも順に並べて添字をつけたものを $\{{\bf x}_1, {\bf x}_2\, \cdots , {\bf x}_{n}\}$ とする。ここで、以下のような写像 $\phi : V \to \mathbb{R}^n$ を定義する。

$\displaystyle{ \phi({\bf x}_i) = {\bf e}_i\ (i=1, 2, \cdots , n) }$

ここで、 ${\bf e}_i$ は $\mathbb{R}^n$ の標準基底である。これによって $T$ の表現行列 $A$ が定まる。 $A$ によって実現される $\mathbb{R}^n$ 上の線形変換を $T_A$ とすると、以下のような可換図式が得られる。

$\displaystyle{ \xymatrix{ V \ar[r]^{T} \ar[d]^{\phi} & V \ar[d]^{\phi} \\ \mathbb{R} \ar[r]^{T_A} & \mathbb{R} } }$

これを成立させるためには、 $A$ はどのような行列であれば良いだろうか？これを一般的な状況で説明するのは非常に煩雑なので、ここでは以下のような具体的な設定の元で議論を進める。

$n = 6$
固有値は $\lambda_1, \lambda_2$ の2つであり、それぞれの代数的重複度は5, 1である。
について、以下の2つのJordan鎖がある。
- - $(T - \lambda_1 I){\bf x}_1 = 0$
  - $(T - \lambda_1 I){\bf x}_2 = {\bf x}_1$
  - $(T - \lambda_1 I){\bf x}_3 = {\bf x}_2$
- - $(T - \lambda_1 I) {\bf x}_4 = 0$
  - $(T - \lambda_1 I){\bf x}_5 = {\bf x}_4$
に対応する固有値ベクトルはとする。
- $(T - \lambda_2 I) {\bf x}_6 = 0$

まず、上で示した可換図式により $T_A = \phi \circ T \circ \phi^{-1}$ となる。また、具体的な設定の中で示した式を変形すると以下のようになる。

$\displaystyle{ \begin{eqnarray} &&T{\bf x}_1 = \lambda_1 {\bf x}_1 \\ &&T{\bf x}_2 = {\bf x}_1 + \lambda_1 {\bf x}_2 \\ &&T{\bf x}_3 = {\bf x}_2 + \lambda_1 {\bf x}_3 \\ &&T{\bf x}_4 = \lambda_1 {\bf x}_4 \\ &&T{\bf x}_5 = {\bf x}_4 + \lambda_1 {\bf x}_5 \\ &&T{\bf x}_6 = \lambda_2 {\bf x}_6 \end{eqnarray} }$

さらに、これらに $\phi$ を作用させると以下のようになる。

$\displaystyle{ \begin{eqnarray} &&(\phi \circ T){\bf x}_1 = \lambda_1 {\bf e}_1 \\ &&(\phi \circ T){\bf x}_2 = {\bf e}_1 + \lambda_1 {\bf e}_2 \\ &&(\phi \circ T){\bf x}_3 = {\bf e}_2 + \lambda_1 {\bf e}_3 \\ &&(\phi \circ T){\bf x}_4 = \lambda_1 {\bf e}_4 \\ &&(\phi \circ T){\bf x}_5 = {\bf e}_4 + \lambda_1 {\bf e}_5 \\ &&(\phi \circ T){\bf x}_6 = \lambda_2 {\bf e}_6 \end{eqnarray} }$

上式に ${\bf x}_i = \phi^{-1}({\bf e}_i)$ を代入することで、 $T_A$ は以下のような変換であることが分かる。

$\displaystyle{ \begin{eqnarray} &&T_A{\bf e}_1 = \lambda_1 {\bf e}_1 \\ &&T_A{\bf e}_2 = {\bf e}_1 + \lambda_1 {\bf e}_2 \\ &&T_A{\bf e}_3 = {\bf e}_2 + \lambda_1 {\bf e}_3 \\ &&T_A{\bf e}_4 = \lambda_1 {\bf e}_4 \\ &&T_A{\bf e}_5 = {\bf e}_4 + \lambda_1 {\bf e}_5 \\ &&T_A{\bf e}_6 = \lambda_2 {\bf e}_6 \end{eqnarray} }$

よって、 $A$ は以下のような行列であることが分かる。

$\displaystyle{ A = \left( \begin{array}{cccccc} \lambda_1 & 1 & 0 & 0 & 0 & 0 \\ 0 & \lambda_1 & 1 & 0 & 0 & 0 \\ 0 & 0 & \lambda_1 & 0 & 0 & 0 \\ 0 & 0 & 0 & \lambda_1 & 1 & 0 \\ 0 & 0 & 0 & 0 & \lambda_1 & 0 \\ 0 & 0 & 0 & 0 & 0 & \lambda_2 \end{array} \right) }$

これはJordan標準形そのものである。結局、線形変換の基底として広義固有ベクトルを選んだ場合の表現行列こそが、Jordan標準形の正体なのである。

広義固有ベクトルへの座標変換

最後に、適当な行列をJordan標準形に変形するとはどういうことなのか、その意味を考えてみよう。

Jordan標準形になっていない $n$ 次正方行列 $B$ を考える。 $B$ をあるベクトル空間 $V$ 上の線形変換 $T$ の表現行列であると考えると、これは基底として広義固有ベクトル以外のものを選んだ場合であると解釈できる。これをJordan標準形に変形する事は、基底を広義固有ベクトルに取り替えることを意味する。

$B$ が表現行列となるような $V$ の基底を ${\bf y}_1, {\bf y}_2, \cdots, {\bf y}_n$ とし、写像 $\psi: V \to \mathbb{R}^n$ を以下のように定める。

$\displaystyle{ \psi({\bf y}_i) = {\bf e}_i\ (i=1, 2, \cdots , n) }$

すると、以下の可換図式が得られる。

$\displaystyle{ \xymatrix{ \mathbb{R} \ar[r]^{T_B} & \mathbb{R} \\ V \ar[r]^{T} \ar[d]^{\phi} \ar[u]_{\psi} & V \ar[d]^{\phi} \ar[u]_{\psi} \\ \mathbb{R} \ar[r]^{T_A} & \mathbb{R} } }$

これより、広義固有ベクトルを基底とした場合の変換 $T_A$ は以下のようにして得られる。

$\displaystyle{ \begin{eqnarray} T_A &=& \phi \circ \psi^{-1} \circ T_B \circ \psi \circ \phi^{-1} \\ &=& (\psi \circ \phi^{-1})^{-1} \circ T_B \circ (\psi \circ \phi^{-1}) \end{eqnarray} }$

$\psi \circ \phi^{-1}$ は基底の取り替えを表す写像である。この写像の表現行列を $P$ とすると、よく見慣れたJordan標準形への変換式 $A = P^{-1}BP$ が得られる。

ちなみに、 $P$ を具体的に求めると、これは広義固有ベクトルを列ベクトルとして順に並べた行列になっている。そのようになる理由は説明すると長くなるので、私の体力の都合により割愛する。

まとめ

本稿ではJordan標準形の理論的側面について述べた。結論として、Jordan標準形とは広義固有ベクトルを基底としたときの線形変換の表現行列であることが分かった。また、Jordan標準形への変換とは、基底を広義固有ベクトルに取り替える操作であることが分かった。

本稿ではJordan標準形の応用面について触れることが出来なかった。これについてはまた機会があれば調べてみたい。

参考

[1] ときわ台学/固有値論/一般固有値問題，一般固有空間，ジョルダン標準形
[2]

線型代数入門 (基礎数学1)

作者: 齋藤正彦
出版社/メーカー: 東京大学出版会
発売日: 1966/03/31
メディア: 単行本
購入: 4人クリック: 102回
この商品を含むブログ (47件) を見る

[3] 線形代数Ｉ/広義固有空間の構造とジョルダン標準形 - 武内＠筑波大

2018-09-17

固有ベクトルの観点から線形変換を掌握する

線形代数

線形代数において、固有値と固有ベクトルの話題は花形である。これらは理論的に美しいだけでなく、応用上様々な場面に登場し、重宝されている。中でも固有値が力を発揮するのは、やはり行列の対角化を行う時であろう。すなわち、n次正方行列 $A$ がある条件を満たすと、ある変換行列 $P$ が存在して $P^{-1} A P$ が対角行列になるのである。

では、対角化が可能となるための「ある条件」とは一体なんだろう？それは、固有値の (固有方程式の解としての) 重複度と、その固有値に対応する固有空間の次元が全ての固有値に対して等しいことである。これは対角化可能であるための必要十分条件となっている。

そう、固有値というのは実に厄介で、例えば固有方程式を解いた結果、 $\lambda$ が重解として得られたとしても、 $\lambda$ に対応する固有空間は必ずしも2次元になるとは限らないのである。

以上の事実から、行列は固有値の数と重複度、および固有ベクトルの次元によって、いくつかの種類に分類出来そうな気がしてくる。そのように分類された行列のクラスは、それぞれ何か特徴的な性質を持つのだろうか？

本稿ではこのぼんやりとした疑問の答えに迫るべく、前半は対角化について調べ、後半は具体例を頼りに固有値・固有空間が反映する行列の性質について調べてみようと思う。

対角化とは何か？

まずは対角化について考えてみよう。対角化と言えば $P^{-1} A P$ という変換の仕方が特徴的である。これの意味するところを考えることによって、対角化について理解を深めていこう。

線形変換の表現

行列は線形変換と密接な関係にある。すなわち、任意の線形変換は行列によって表現することができる。 $P^{-1} A P$ の意味を理解する上では、Aをある線形変換の表現行列と捉えるのがよい。

まず、線形変換の定義を[1]より引用する。

線形変換

$K$ 上の線形空間 $V$ から $K$ 上の線形空間 $V'$ への写像 $T$ が二条件
$\displaystyle{ \begin{eqnarray} && T({\bf x} + {\bf y}) = T({\bf x}) + T({\bf y}) \\ && T(a{\bf x}) = aT({\bf x}) \end{eqnarray} }$
を充すとき, $T$ を $V$ から $V'$ への線形写像と言う.
・・・ (中略) ・・・
とくに, $V$ から $V$ 自身への線形写像を, $V$ の線形変換と言う.

$V$ の次元をnとする。 $V$ の1つの基底を ${\bf e}_1, {\bf e}_2, \cdots , {\bf e}_n$ とすると、 $V$ の任意の元は以下のように基底の線型結合で表すことができる。

$\displaystyle{ {\bf v} = \sum_{i=1}^n x_i {\bf e}_i }$

ここで、 ${\bf v}$ を $V$ の線形変換 $T$ で写したものはやはり $V$ の元となる。よって以下のように同じ基底で表すことができる。

$\displaystyle{ T({\bf v}) = \sum_{i=1}^n x_i' {\bf e}_i }$

すると、変換前後の係数の間には以下のような関係が成立する。

$\displaystyle{ \left( \begin{array}{c} x_1' \\ x_2' \\ \vdots \\ x_n' \end{array} \right)_e = \left( \begin{array}{cccc} a_{11} & a_{12} & \ldots & a_{1n} \\ a_{21} & a_{22} & \ldots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \ldots & a_{nn} \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right)_e }$

上の行列 $A = (a_{ij})$ が線形変換 $T$ の表現行列である。両辺のベクトルの右下にeと書いたのは、これが基底 ${\bf e}_i$ によるベクトルであることを示している。表現行列を簡潔に $A$ と表した式は以下のようになる。

$\displaystyle{ \left( \begin{array}{c} x_1' \\ x_2' \\ \vdots \\ x_n' \end{array} \right)_e = A \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right)_e }$

基底の取り替え

次に、線形変換 $T$ を $V$ の別の基底で表現するとどうなるか考えてみよう。そのような基底を ${\bf f}_1, {\bf f}_2, \cdots , {\bf f}_n$ とする。すると、先ほどと同じように変換前後のベクトルは以下のように表される。

$\displaystyle{ {\bf v} = \sum_{i=1}^n y_i {\bf f}_i }$

$\displaystyle{ T({\bf v}) = \sum_{i=1}^n y_i' {\bf f}_i }$

基底の取り替えにおいて重要な考え方は、元の基底に対する表現行列 $A$ をできる限り活かすことである。そのために、 $A$ ではなく入力として与えるベクトルの方を変換することを考える。これを実現するためには、以下のような行列 $P$ を利用する。

$\displaystyle{ \left( \begin{array}{c} x_1 \\ x_2 \\ \vdots \\ x_n \end{array} \right)_e = P \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right)_f }$

このような行列 $P$ を基底の取り替え行列と呼ぶ。これより、 $T$ は以下のように表現できる。

$\displaystyle{ P \left( \begin{array}{c} y_1' \\ y_2' \\ \vdots \\ y_n' \end{array} \right)_f = A P \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right)_f }$

異世界を旅して帰還する

基底の変換後の式の両辺に $P^{-1}$ をかけると以下のようになる。

$\displaystyle{ \left( \begin{array}{c} y_1' \\ y_2' \\ \vdots \\ y_n' \end{array} \right)_f = P^{-1}A P \left( \begin{array}{c} y_1 \\ y_2 \\ \vdots \\ y_n \end{array} \right)_f }$

$P^{-1}A P$ が基底 ${\bf f}_i$ を用いた場合の $T$ の表現行列となる。

普通はここで話はおしまいである。基底を変換した後の表現行列が求められたのだから、めでたしめでたしというわけだ。しかし、 $P^{-1}A P$ という形の変換は線形代数に限らずいろんな数学の分野でよく見かける形であり、もう一段抽象的なレベルでの意味がある。これを少し掘り下げて考えてみよう。

まず、行列 $A$ は基底 ${\bf e}_i$ の世界で表現されるベクトルに対して適用可能な行列である。これを別の世界のベクトル、すなわち、基底 ${\bf f}_i$ で表現されるベクトルに適用したいと思ったら、一度 ${\bf f}_i$ の世界から ${\bf e}_i$ の世界に移動する必要がある。すると、移動後のベクトルには $A$ を適用することができる。最後に、 $A$ を適用して得られたベクトルを元の世界に戻してやれば、結局 ${\bf f}_i$ の世界のベクトルに $A$ を適用出来たことになる。以下に図を示す。

f:id:peng225:20180917213651p:plain

このように $P^{-1}, P$ で変換を挟み込むことは、その変換が適用可能な世界へと移動し、変換が終わったら元の世界に戻してやるような効果がある。これが線形変換の表現行列を別の基底で表したベクトルに適用できるからくりである。

対角化は座標変換

ここまで座標変換の話ばかりしてきた訳だが、座標変換と対角化の関係について述べておく。ずばり対角化とは、線形変換の基底を固有ベクトルから成る基底に変換する操作のことである。そして、対角化可能であるとは、基底となるような固有ベクトルの組が存在することを意味している。

逆に、もし固有値の重複度より固有ベクトルの次元が小さいものが存在すると、対象となる線形空間の次元に対して固有ベクトルの本数が足りなくなってしまい、基底を作れなくなる。この場合は対角化不可能となる。

基底の取り替え行列 $P$ としては、 $A$ のn個の列固有ベクトル ${\bf p}_1, {\bf p}_2, \cdots, {\bf p}_n$ を並べた行列を取れば対角化できる。計算の過程を以下に示す。

$\displaystyle{ \begin{eqnarray} P^{-1} A \left( \begin{array}{cccc} {\bf p}_1 & {\bf p}_2 & \cdots & {\bf p}_n \end{array} \right) &=& P^{-1} \left( \begin{array}{cccc} \lambda_1 {\bf p}_1 & \lambda_2 {\bf p}_2 & \cdots & \lambda_n {\bf p}_n \end{array} \right) \\ &=& \left( \begin{array}{cccc} \lambda_1 & 0 & \ldots & 0 \\ 0 & \lambda_2 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \lambda_n \end{array} \right) \end{eqnarray} }$

固有値・固有ベクトルは何に紐づくものか？

ここで注意点を1つ。世の中ではよく「行列の固有値」とか「行列の固有ベクトル」という言い方をする。これは別に間違っていない。実は、線形変換に対しても固有値・固有ベクトルという概念が存在する。

まず固有値についてだが、行列の固有値は $P^{-1}, P$ のサンドイッチ変換に対して不変となる。そのため、線形変換の固有値とは「任意の基底に対する表現行列の固有値である」と定義しておけば、それは一意に定まる。

多様体や微分幾何学なんかが代表的であるが、座標系に依存しない量や概念というのは、その数学的対象の本質的な性質を表していると考えられるため、一般にとても重要である。

一方、固有ベクトルは座標系の取り方によって変わる。このように、固有値と固有ベクトルでは不変となる範囲が異なるので注意が必要である。

固有値・固有空間に着目した線形変換の分類

ここまでで固有値・固有ベクトルとかなりお友達になれたはずなので、本題に入ろう。線形変換が固有値の数や重複度、また固有空間の次元とどのような関係にあるのかを調べてみよう。簡単のために、話を $\mathbb{R}$ 上2次の正方行列に絞る。

$\mathbb{R}$ 上2次の正方行列は以下のように分類できる。

異なる2つの実固有値を持つ場合
1つの実固有値を持ち、固有空間の次元が2の場合
1つの実固有値を持ち、固有空間の次元が1の場合
実固有値を持たない場合

以下では分類された各項目の線形変換がどのような振る舞いをするのか、例を用いながら調べてみる。

異なる2つの実固有値を持つ場合

例として以下の行列を考える。

$\displaystyle{ A = \left( \begin{array}{cc} 1 & -1/4 \\ -1/2 & 5/4 \\ \end{array} \right) }$

$A$ の固有値・固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[1, -1/4], [-1/2, 5/4]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(3/2, [   # 固有値1つ目
(1, -2)   # 固有ベクトル
], 1)     # 重複度
(3/4, [   # 固有値2つ目
(1, 1)    # 固有ベクトル
], 1)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002223p:plain

ここで、黒の矢印は基底、グレーの点線矢印は基底に $A$ を作用させて得られるベクトルである。赤の矢印は固有ベクトル、オレンジの点線矢印は固有ベクトルに $A$ を作用させて得られるベクトルである。また、水色の矢印は、矢印の根元の位置ベクトルが、 $A$ によって矢印の先端に写ることを意味している。

これを見ると、Aの作用によって空間が2つの固有ベクトルの方向に伸び縮みしているのが分かる。これより、このタイプの線形変換は固有ベクトル方向の拡大・縮小を表すと考えられる。ただし、伸びるのか縮むのか、またその程度がどれくらいであるかは固有値によって決まる。

1つの実固有値を持ち、固有空間の次元が2の場合

例として以下の行列を考える。

$\displaystyle{ A = \left( \begin{array}{cc} 2 & 0 \\ 0 & 2 \\ \end{array} \right) }$

$A$ の固有値・固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[2, 0], [0, 2]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(2, [     # 固有値
(1, 0),   # 固有ベクトル1つ目
(0, 1)    # 固有ベクトル2つ目
], 2)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002458p:plain

これを見ると、Aの作用によって空間が原点を中心に拡大されているのが分かる。これより、このタイプの線形変換は拡大・縮小を行うものと考えられる。

一点補足だが、このケースの線形変換の表現行列は単位行列の定数倍しかあり得ない。 $A$ のただ1つの固有値を $\lambda$ とすると、 $A$ を対角化した行列は $\lambda E$ となる。対角化のための変換行列を $P$ とすると、結局以下のようになる。

$\displaystyle{ \begin{eqnarray} P^{-1} A P &=& \lambda E \\ A &=& P (\lambda E) P^{-1} \\ A &=& \lambda E \end{eqnarray} }$

よって $A$ は単位行列の固有値倍となる。

1つの実固有値を持ち、固有空間の次元が1の場合

例として以下の行列を考える。

$\displaystyle{ A = \left( \begin{array}{cc} 11/6 & -1/3 \\ 1/3 & 7/6 \end{array} \right) }$

$A$ の固有値・固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[11/6, -1/3], [1/3, 7/6]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(3/2, [   # 固有値
(1, 1)    # 固有ベクトル
], 2)     # 重複度

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002520p:plain

これを見ると、原点を通る固有ベクトルによって引かれる直線より下側では右向きの、上側では左向きの流れがあるように見える。これより、このタイプの線形変換は固有ベクトルと平行な方向に歪みを与える働きがあると考えられる。

実固有値を持たない場合

例として以下の行列を考える。

$\displaystyle{ A = \left( \begin{array}{cc} 1 & -1/2 \\ 1/2 & 1 \end{array} \right) }$

$A$ の固有値・固有ベクトルは以下のようになる。

sage: A = matrix(QQ, [[1, -1/2], [1/2, 1]])
sage: for var in A.eigenvectors_right():
....:     print var
....:     
(1 - 0.50000000000000000?*I, [(1, 1*I)], 1)    # 固有値が虚数
(1 + 0.50000000000000000?*I, [(1, -1*I)], 1)

この線形変換を可視化したものを以下に示す。

f:id:peng225:20180913002536p:plain

これを見ると、螺旋のような渦巻きが現れている。これより、このタイプの線形変換は回転と拡大・縮小の組み合わせになっていると考えられる。

まとめ

以上、線形変換についていろいろと調べてみた。その結果、線形変換は固有値と固有ベクトルによっていくつかの種類に分類できそうだということが分かった。本稿では2次正方行列で表現される線形変換しか詳細に調べなかったが、3次以上についても類似の考察が可能だろうと推測される。3次までなら頑張って可視化できないこともないかもしれないので、機会があればトライしてみたい。

参考文系

[1]

線型代数入門 (基礎数学)

作者:齋藤正彦
出版社/メーカー: 東京大学出版会
発売日: 1966/03/31
メディア: 単行本

2018-07-31

区間推定の理論と実践

統計学

統計学における推定の考え方は仮説検定と双璧を成す重要な理論であり、初学者がまず目指すべき最高到達点である。推定は強力な手法であるが故に、恐らく大抵の教科書では記載があるだろうし、ググればいくらでもWebページがヒットする。

そんな状況下ではあるが、推定の中でも特に区間推定は個人的に興味をそそられる内容であり、ぜひともブログを書き起こしてみたくなった。そこで、本稿ではできるだけ二番煎じにならないよう、私が得た区間推定に対する理解について説明してみたいと思う。

本稿では正規母集団から抽出した標本データを用いた母平均の区間推定のみをターゲットとする。また、適当なサンプルデータを用いて実際に区間推定を行うことで、区間推定の理論に対して確かな実感を得てみたいと思う。

なお、本稿を執筆するに辺り、全体的に[1][2]を参考にした。

母平均の区間推定

母分散が既知の場合

正規母集団 $N(\mu,\ \sigma^2)$ に従う独立な確率変数 $X_1,\ X_2,\ \cdots .\ X_n$ の標本平均 $\overline{X}$ を考える。このとき、 $\overline{X}$ 自体も確率変数となっており、これが従う確率分布は $N(\mu,\ \sigma^2/n)$ となる。また、これを標準化した確率変数 $Z = \frac{\overline{X} - \mu}{\sigma/\sqrt{n}}$ は、当然だが標準正規分布 $N(0,\ 1)$ に従う。

今、母分散 $\sigma^2$ は分かっているが、母平均 $\mu$ の値が不明な状況を考える。もちろん、実際にはこんなことはまずあり得ないのだが、このケースを考えることは後々の議論において意味がある。

この時、 $\mu = E(\overline{X})$ という関係が成り立つので、標本平均をそのまま母分散だと考えてもあながち間違いではない。しかし、当然だが $\overline{X}$ の値にはぶれがあり、どんなサンプルを抽出するかによって値が変わってしまう。しかも、極端なケースでは $\overline{X}$ はとんでもなく大きくなったり小さくなったりすることがあり得る。そのため、 $\overline{X}$ の値だけをみても、母平均について何も確実なことは言えない。

では、何か確率的な議論はできないだろうか？例えば「95%の確率で $\mu$ はAからBの間に入ります」というようなことが言えないだろうか？

結論から言うとこれは可能である。標準正規分布において、それより大きい値を取る確率が $\alpha\ (0 < \alpha < 1)$ となるような点 ( $100 \alpha$ パーセント点) を $Z_\alpha$ と表すと、先ほど導入した確率変数Zは確率 $1 - \alpha$ で以下の不等式を満たす。

$\displaystyle{ -Z_{\alpha/2} < Z < Z_{\alpha/2} }$

2点注意事項がある。まず、これは正規分布が左右対称だから言えることである。さもなければ、 $-Z_{\alpha/2}$ の部分はそれより小さい値を取る確率が $\alpha/2$ となるような点を別途求めて、それと置き換えてやる必要がある。

次に、 $\alpha$ を2で割っている理由だが、これは上側の $\alpha/2$ と下側の $\alpha/2$ を両方排除し、残りが $1 - \alpha$ となるようにするためである。

この不等式を変形すれば、以下のように $\mu$ が確率 $1 - \alpha$ で取る値の範囲が分かる。

$\displaystyle{ \overline{X} - \frac{Z_{\alpha/2} \sigma}{\sqrt{n}} < \mu < \overline{X} + \frac{Z_{\alpha/2}\sigma}{\sqrt{n}} }$

母分散が未知の場合

先ほども述べた通り、母平均は未知なのに母分散は分かっているということは普通はあり得ない。実際に得られるのは不偏分散 $s^2$ だけというケースがほとんどである。不偏分散の定義式を以下に示す。

$\displaystyle{ s^2 = \frac{1}{n-1} \sum_{i=1}^n (X_i - \overline{X})^2 }$

不偏分散に対しては $\sigma^2 = E(s^2)$ という関係式が成り立つので、これを母分散の代わりに用いることはそれほど不自然ではない。しかし、やはりサンプル数が少ない場合には母分散とかけ離れた値を取ることがあり得る。そのため、先ほど考えた確率変数Zの式に対して単に $\sigma = s$ を代入すると、これはもはや正規分布には従わない。

では、これはどのような分布になるのだろうか？それを調べるために、Zの式に $\sigma = s$ を代入した式を以下のように変形してみよう。

$\displaystyle{ \begin{eqnarray} t &=& \frac{\overline{X} - \mu}{\frac{s}{\sqrt{n}}} \\ &=& \frac{\overline{X} - \mu}{\frac{\sigma}{\sqrt{n}} \cdot \frac{s}{\sigma}} \\ &=& \frac{Z}{\frac{\sqrt{\frac{1}{n-1} \sum_{i=1}^n (X_i - \overline{X})^2}}{\sigma}} \\ &=& \frac{Z}{\sqrt{\frac{1}{n-1} \sum_{i=1}^n \left(\frac{X_i - \overline{X}}{\sigma}\right)^2}} \end{eqnarray} }$

ここで、分母に $\sum_{i=1}^n \left(\frac{X_i - \overline{X}}{\sigma}\right)^2$ という式が現れた。この式自体も確率変数であり、これは標準正規分布に従う確率変数の二乗和を表す。このような確率変数が従う確率分布は自由度n-1の $\chi^2$ 分布と呼ばれる。

まとめると、標準正規分布に従う確率変数Zと自由度n-1の $\chi^2$ 分布に従う確率変数Wによって、元々の確率変数tは以下のように表される。

$\displaystyle{ t = \frac{Z}{\sqrt{\frac{W}{n-1}}} }$

このような確率変数tが従う確率分布は自由度n-1のt分布と呼ばれる。結論として、母分散未知のケースにおける母平均の推定にはt分布を使用する必要があることが分かった。

自由度n-1のt分布を $t(n-1)$ と書くと、母分散既知のケースと同様の計算を行うことで、 $\mu$ は確率 $1 - \alpha$ で以下の不等式を満たすことが分かる。

$\displaystyle{ \overline{X} - \frac{t_{\alpha/2}(n-1) s}{\sqrt{n}} < \mu < \overline{X} + \frac{t_{\alpha/2}(n-1) s}{\sqrt{n}} }$

ただし、 $t_{\alpha/2}(n-1)$ は自由度n-1のt分布の $100 \alpha/2$ パーセント点を表す。ここでもt分布が左右対称であることを利用している。

サンプル数に関する注意

t分布は自由度が大きくなると正規分布に収束していく。そのため、t分布を用いる必要があるのはサンプル数が少ない場合に限る。実用上はサンプル数が30以上あれば正規分布を使用しても差し支えないようである[1]。

区間推定の実践

ここまでに述べたことを実践してみよう。統計処理をするにはpythonかRが簡単そうであるが、Android上で環境を構築できたのでpythonを使うことにする[3]。また、データセットとしてIrisは見飽きたので、ここではAbarone (アワビ) を使用している[4]。

Abaroneデータセットには、1列目に性別、3列目に直径の情報が入っている。今回はオスのアワビの直径データについて、母平均の区間推定を行ってみた。サンプル数は敢えて少なめに10とし、t分布を使った場合と正規分布を使った場合で比較を行った。もちろん、正しいのはt分布を使う方である。

以下に今回作成したソースコードを示す。作成の際には[5]を参考にさせて頂いた。日本語コメントは本稿のために後から付したものであり、このままでは実行できないので注意されたい。

import pandas as pd
from scipy import stats
import math

# 全部の列に名前を付けるのは諦めた...
abalone = pd.read_csv('file:///storage/emulated/0/Download/abalone.data', 
                 header=None, names=['sex', 'length', 'diameter', 'height', 'a', 'b', 'c', 'd', 'e'])

# オスの情報だけ抽出
male = abalone[abalone['sex'] == 'M']

n = 10
# 自由度n-1のt分布の95%信頼区間を取得
tleft, tright = stats.t.interval(alpha=0.95, df=n-1)
# 正規分布の95%信頼区間を取得
nleft, nright = stats.norm.interval(alpha=0.95)

mean = male.iloc[:n]['diameter'].mean()
# デフォルトで不偏分散の平方根が計算される
stddev = male.iloc[:n]['diameter'].std()

# 本文中の式と若干符号が異なるので注意
tleft = mean + tleft*stddev/math.sqrt(n)
tright = mean + tright*stddev/math.sqrt(n)

nleft = mean + nleft*stddev/math.sqrt(n)
nright = mean + nright*stddev/math.sqrt(n)

print("mean = {:.3f}, stddev = {:.3f}".format(mean, stddev))
print("t-dist: left = {:.3f}, right = {:.3f}".format(tleft, tright))
print("n-dist: left = {:.3f}, right = {:.3f}".format(nleft, nright))

実行結果を以下に示す。

mean = 0.339, stddev = 0.046
t-dist: left = 0.306, right = 0.372
n-dist: left = 0.310, right = 0.368

実行結果から分かるように、正規分布を使った場合は区間が狭めに出力されてしまう。そのため、正規分布での結果を元に考察を行うと判断を誤る恐れがある。

最後に特大の注意事項がある。今回はオスのアワビの直径データが正規分布に従うと仮定して分析を行ったが、本来であればこの仮定の妥当性は良く確かめる必要がある。誤った仮定を元に一生懸命計算したとしても、その結果には何の意味もないどころか、誤った帰結へと導かれる危険性がある。

データが正規分布に良く当てはまっているかどうかを判断する方法について私には知見がない。軽くググってみたらいろいろと情報がありそうだったので、興味のある方は調べてみると面白いかもしれない。

その他の話題

本稿では正規母集団から抽出したサンプルから母平均の区間推定を行った。この他にも、例えば母分散の区間推定を行う方法などは有名である[1]。また、実用上のタスクとしては正規母集団以外の母集団を考えることもあるだろう。しかし、どんな場合であれ、基本的な区間推定の考え方は本稿に述べたものと大きく異なることはないであろう。

まとめ

以上、区間推定に関する理論と実践について述べた。統計学は応用先が非常にたくさんある分野なので、理解を深めておいて損はない。また機会があれば、多変量解析などの進んだ内容も勉強し直したい。

参考

[1]

統計学入門 (基礎統計学?)

作者: 東京大学教養学部統計学教室
出版社/メーカー: 東京大学出版会
発売日: 1991/07/09
メディア: 単行本
購入: 158人クリック: 3,604回
この商品を含むブログ (79件) を見る

[2] 統計学の時間 | 統計WEB
[3] https://play.google.com/store/apps/details?id=ru.iiec.pydroid3&hl=ja&referrer=utm_source%3Dgoogle%26utm_medium%3Dorganic%26utm_term%3Dpydroid+3&pcampaignid=APPU_1_4VFeW8v9NovS8wXKsq_4Cw
[4] UCI Machine Learning Repository: Abalone Data Set
[5] Python: SciPyを使った統計的推定 - け日記

2018-07-22

確率分布族の再生性まとめ

確率分布はたくさんある

統計学を勉強していると、世の中には本当にたくさんの確率分布が存在する事に気付く。現実世界の現象に統計学を応用しようと考えるとき、その現象にはどのような確率分布が良く当てはまるのかを考えたくなる時がある。そのためには、数ある確率分布がそれぞれどういう現象をうまく表現できるかということを知っておく必要がある。

そこで、本稿では基本的な確率分布の特徴についてまとめてみ・・・ようと思ったのだが、途中までやりかかったところですでに素晴らしいまとめ[1]が存在することに気づいた。

劣化コピーを作っても意味がないので、ここでは[1]に記載が少ない再生性に関する事項についてのみまとめてみる。ただし、さすがに[1]に記載されている全ての分布について調べるのはしんどいので、ここでは私が個人的に興味を持った分布についてのみ紹介する。詳細は[2][3][4][5][6]等を参照のこと。

離散型

名称	再生性	$X_i\ (i=1,\ 2)$	$X_1 + X_2$
二項分布	有り	$Bi(n_i,\ p)$	$Bi(n_1 + n_2,\ p)$
ポアソン分布	有り	$Po(\lambda_i)$	$Po(\lambda_1 + \lambda_2)$

連続型

名称	再生性	$X_i\ (i=1,\ 2)$	$X_1 + X_2$
正規分布	有り	$N(\mu_i,\ \sigma_i^2)$	$N(\mu_1 + \mu_2,\ \sigma_1^2 + \sigma_2^2)$
対数正規分布	有り	$\Lambda(\mu_i,\ \sigma_i^2)$	$\Lambda(\mu_1 + \mu_2,\ \sigma_1^2 + \sigma_2^2)$
指数分布	無し	-	-
ガンマ分布	有り	$Ga(\alpha_i,\ \lambda)$	$Ga(\alpha_1 + \alpha_2,\ \lambda)$
コーシー分布	有り	$Ca(\alpha_i,\ \beta_i)$	$Ca(\alpha_1+\alpha_2,\ \beta_1+\beta_2)$
$\chi^2$ 分布	有り	$\chi^2(m_1)$	$\chi^2(m_1 + m_2)$
t分布	無し？	-	-