気になる確率の話

火星に猫はいるか？

Feb 09, 2012

火星に生物はいるだろうか？生物のいる確率はどのくらいだろうか？順を追って求めてみよう。

まずは火星に人間がいる確率を $P_\text{human}$ としよう。例えば… $P_\text{human}=0.0000001\%$ くらいかな。

よし、次に火星に犬がいる確率を $P_\text{dog}$ としよう。なんとなく犬は人間よりいそうなんで、 $P_\text{dog}=0.000001\%$ とでもしてみるか。

どんどん続けよう。次に火星にアメンボがいる確率を $P_\text{water strider}$ としよう。これはきっと人間なんかよりずーっといる確率が高いだろう。よし $P_\text{water strider}=0.001\%$ だ。

ここで中間報告。では、火星に人間も犬もアメンボもいない確率はいくつになるだろうか。それは、”（人間がいない）かつ（犬がいない）かつ（アメンボがいない）”確率となるので、次式で与えられる。

$P_\text{no life}=(1-P_\text{human})(1-P_\text{dog})(1-P_\text{water strider})$ $=99.9999999\%\times 99.999999\%\times 99.999\%=99.9989989\%$

なるほど…しかし、まだ他にも調査すべき生物はたくさんある！象がいる確率、マントヒヒがいる確率、ムカデがいる確率、カピバラがいる確率…どれもきっと確率は0ではないだろう。

ある生物がいる確率を $P_{i}$ とする。 $i$ は生物の種別を表していている。ものの本によると生物は1億種という説もある。それらが火星にいる確率が平均して $0.000001\%$ とすると、火星にまったく生物がいない確率は、

$P_\text{no life}=\Pi_{i=1}^{100000000}(1-P_{i})=(99.999999\%)^{100000000}\sim 37\%$

となる。ちょっと待てよ、、、ということは、火星に何か生物がいる確率は、

$1-P_\text{no life}\sim 63\%$

これは驚いた！火星にはきっと何かいるぞ！

上記の推論はどこがおかしいか分かりますが？

男の子？女の子？

Dec 20, 2011

「二人の子供問題」という有名な確率のトピックです。

昼休み。食事も終えて、職場の仲間とまったりとコーヒーを飲んでいます。誰ともなしに家族の話題になりました。

独りモノのあなたな職場の仲間Aさんにこんな質問をしました。

質問1：子供は何人いますか？
回答1：二人です

では、ここで問題

問1：Aさんの子供が、二人とも男の子の確率は？

あなたはAさんに更に質問をぶつけました。

質問2：では、男の子はいますか？
回答2：はい

なるほど。では、ここで問題

問2：Aさんの子供が、二人とも男の子の確率は？

もう一人Bさんが隣に座っています。Bさんは今日はたまたま子供を連れてきています。男の子です。

あなたはBさんにも質問をしました。

質問3：子供は何人いますか？
回答3：二人です

ふーむ、、、ここで問題

問3：Bさんの子供が、二人とも男の子の確率は？

問2と問3の違いが分かりますか？AさんもBさんも二人子供がいて、少なくとも一人は男の子であることが分かっています。

ちなみに解答は（男の子が生まれる確率と女の子が生まれる確率は等しいとしましょう）

解1：1/4
解2：1/3
解3：1/2

です。

メルセンヌ素数

Dec 20, 2011

確率のネタではないですが、小話をひとつ

NtRandでは”Mersenne-Twister”と呼ばれる一様乱数生成アルゴリズムが採用されています。このアルゴリズムで生成される乱数は $2^{19937}-1$ 、実に6002桁という途方もない長い周期をもっています。

さて、この数は”メルセンヌ素数”と呼ばれる一群の数字のうちの一つで、これが”Mersenne-Twister”の語源になっています。

メルセンヌ素数とは、素数のうちで $2^{n}-1$ という形で書けるもののことです。つまりMarsenne-Twisterの周期は $n=19937$ のメルセンヌ素数です。 2001年5月現在、メルセンヌ素数は47個見つかっています。 $n=19937$ のメルセンヌ素数は24番目のものです。

1644年、マラン・メルセンヌさんが $n=2,3,5,7,13,17,19,67,127,257$ がメルセンヌ素数であると主張しました。残念なことに彼は $n=61,89,107$ を見逃していました。

更に残念なことに、、、

1903年、アメリカで行われた数学会でフランク・ネルソン・コールなる人物が発表のため登壇しました。タイトルは「大きな数の素因数分解」。

彼はまず黒板に $2^{67}$ を書き下しました。

$2^{67}=147,573,952,589,676,412,928$

そこから1を引き、そして次に

$193,707,721\times761,838,257,287$

の計算を行いました。その結果は、、、 $147,573,952,589,676,412,927$ ! その間彼は一言も口をきかず、静かに席に戻りました。その後会場は万雷の拍手が沸いたそうです。

メルセンヌさんの主張から250年以上経って、 $2^{67}-1$ が素数でないことが分かったのです。

今では $n=127$ も素数ではないことが分かっています。

某国、華やかなカジノの街。
ルーレット、スロットマシーン、カード…有象無象どもが一攫千金を目論んで目を血走らせて必死にもがいている。
さて今日はどのゲームでひと稼ぎしようか…ふと片隅にあるゲームに目がとまる。

「サイコロを転がして出た目が1なら1円、2なら2円…出た目と同じ額を賞金として進呈！」

なんとも幼稚なギャンブルだが、まぁいいや。ちょっとやってみるか。参加費用は4円か。サイコロなんだから5か6を出せばもうかるじゃないか、チョロイチョロイ。
…1時間後…
すっかりスッカラカンになった彼は道端でボロ雑巾のようになっていた。

もちろん、確率を熟知したあなたはこのゲームには参加しないはずですね。そう、期待値の問題です。このゲームの賞金の期待値 $E$ は、

$E=1\times \frac{1}{6}+2\times \frac{1}{6}+3\times \frac{1}{6}+4\times \frac{1}{6}+5\times \frac{1}{6}+6\times \frac{1}{6}=3.5$

です。賞金は 3.5円しか期待できないのに参加費が4円。長い目でみるとどんどんと負けが込んでいくのです（ちなみに1枚300円の宝くじの期待値は大体142円だそうです）。

リベンジだとばかりに再びカジノにやってきた彼。今度は確率の勉強もしっかりしてきた。うーん、どのゲームも期待値を計算するとマイナスのものばかりか（当たり前。ギャンブルは胴元が儲かるようになっている）…なまじ確率の勉強をしたのでどのゲームも割に合わないと感じてしまうなぁ。ふと片隅にあるゲームに目がとまる。
「コイントスゲーム！外れなし！表が出たら賞金は倍！！ただし裏が出たらゲームオーバー」
つまり、

1回目で裏が出たら賞金1円
2回目で裏が出たら賞金2円
3回目で裏が出たら賞金4円
4回目で裏が出たら賞金8円
（以下賞金は倍々）

なるほど面白そうだ。おっといけない、軽々に飛びついては痛い目に逢う。まずは期待値を計算してみよう。
各パターンの確率は、

裏： $1/2$
表裏： $(1/2)^2$
表表裏： $(1/2)^3$
表表表裏： $(1/2)^4$
表表表表裏： $(1/2)^5$
（あとは分かるでしょう）

となるので、期待値は

$E=1\times\frac{1}{2}+2\times\left(\frac{1}{2}\right)^2+4\times\left(\frac{1}{2}\right)^3+8\times\left(\frac{1}{2}\right)^4+\cdots=\frac{1}{2}+\frac{1}{2}+\frac{1}{2}+\cdots=+\infty$

こりゃ凄い、何と期待値は無限大だ！幾らでも儲かるってことか！！で、参加費はいくらだって？1億円か。安い安い！！よ～し、やるぞ！！！…はたして彼の運命やいかに…

どうです？あなたはこのゲームに参加しますか？参加費100円でもやらないのではないでしょうか？（ちなみに100円以上の賞金を手にするには、7回連続して表を出す必要があります）。

でも確かに期待値は無限大…腑に落ちませんね。これが「サンクトペテルブルグのパラドックス」です。

「独立」と「無相関」

May 24, 2010

複数の確率変数があるとき、それらがてんでバラバラ好き勝手に変動しているとき「無相関」という言葉を使います。また「独立」という言葉もあります。
どちらも何となく似たような意味合いなので、混同している（というか気にしていない）場合もあるようですが、もちろんこの2つは異なる定義があります。
では「無相関」と「独立」の関係を見てみましょう。

「独立」なら「無相関」か？

その通り！「独立」の方が強いです。これはまぁ言葉通りのイメージではないでしょうか。「独立」しているのに何らかの関連性があるとは考えにくいです。

「無相関」なら「独立」か？

これは違います！。お互いに関連せずに変動していても「独立」でない場合があります。
これを理解するにはやはり「相関」を正しく理解することが必要です。
相関の定義は、

$\rho_{xy}=\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})(y_{i}-\bar{y})$

です。細かい内容はさておき、重要なのは相関はデータ全体でひとつの値、つまり平均や標準偏差と同じデータの特徴を現す量であることです。
平均は分布の位置、標準偏差は分布の広がり具合です。では相関が何を表す量か？それは下の Flash で確認してください。

相関係数が 0 つまり無相関の場合とそうでない場合のデータプロットの違いが分かりましたか？無相関だと分布が全体としてまん丸になるになりますね。もっとも無相関だからといって分布がまん丸になるとは限りませんが、少なくともまん丸なら無相関です。
では下のようなプロットはどうでしょう？

このプロット (x, y) を生成する確率変数 X と Y は無相関です（まん丸で、どっちに傾けても変化なし）。しかしこの図をみて X と Y が勝手に変動していると認識する人はいないでしょう（たとえば X が 0付近なら、Y は $\pm 1$ 付近しか値がとれない）。

石器時代からあるモンテカルロシミュレーション

May 24, 2010

モンテカルロシミュレーションの第1歩として、100パーセントどの教科書にも載っている例題があります。
どの教科書にも載っていて、このサイトに無いのも悔しいので載せておきます。

一様乱数から円周率を求めよう

1辺が1の正方形を用意します。そこに4分の1円を描いておきます。この正方形内に一様に点を打っていき、（円内の点の数）/（全点の数）を計算すると、これは（4分の1円面積）/（正方形の面積）に近づいていきます。
（4分の1円面積）/（正方形の面積）はすなわち $\pi/4$ なので、（円内の点の数）×4 /（全点の数）は円周率 $\pi$ に近づいていきます。

STDEVP? STDEV?

May 24, 2010

確率や統計にかかわる者が一度は悩む…否、一生悩む問題が「標準偏差を求めるときって、STDEV を使うのか STDEVP を使うのか？」でしょう（STDEV や STDEVP は Excel の関数です。念のため）。この2つの関数の違いは、

STDEVP：（（データを2乗したものの平均）-（データの平均の2乗））/ データ数
STDEV：（（データを2乗したものの平均）-（データの平均の2乗））/ （データ数 – 1）

となっていて、教科書には「データ数から1を引くのは、自由度が1つ減ってどうのこうの…」という説明があり、分かったような分からないような…みなさんはこの2つの標準偏差の違いをしっかり理解していますか？

データの分布：STDEVP

こちらがいわゆる普通の標準偏差です。高校まではこちらの標準偏差しか学ばないでしょう。

あるデータが与えられました。データの羅列だけでは特徴はつかめません。そこでデータを集約して特徴をとらえようと考えます。代表的な特徴値としては、

データの典型的な値は？
データはどれくらいばらついている？

といったところではないでしょうか。データの典型的な値としては平均や中央値、そしてバラツキ度合いが STDEVP で求められる標準偏差です（定義から明らか）。
このように単に与えられたデータの特徴の1つとしての標準偏差が STDEVP です。

データの背後の真の分布：STDEV

1年後の株価を知りたい！これは人間誰しも持つ欲望でしょう。しかし、残念ながら株価は確率的に変動していて正確に知る方法はなさそうです。
株価は確率的に変動する、そして株価の変動（対数収益率）は正規分布に従うというのが一般的な知見です（最近はいろいろと異論があるようですが）。
時々刻々と記録されていく株価の変動は1つの正規分布、つまりある平均 $m$ とある標準偏差 $\sigma$ をもつ正規分布（この2つで正規分布は決定できるというのは常識です）に従っていて、きっと分布の平均 $m$ と標準偏差 $\sigma$ は神様が決めているんでしょう。われわれにはその平均 $m$ と標準偏差 $\sigma$ を正確に知る術はありませんから、なんとかして推定したいと考えます。

例えば株価の変動を毎日毎日記録していきます。1年ごとにデータをまとめて、それぞれのデータセットの平均 $\mu_{i}$ と “標準偏差” $\bar{\sigma}_{i}$ を求めます。ここで “標準偏差” $\bar{\sigma}_{i}$ は 何らかの方法 （具体的には STDEV で求めた標準偏差か STDEVP で求めた標準偏差）に従ったデータのバラツキ具合を表す値です。

結論から言うと、このときの標準偏差は STDEV で求ます。それの根拠は、、、

過去のある1年のデータから”標準偏差” を求めたとします。これは神様が決めた本当の正規分布（無限の過去から無限の未来まですべての株価のデータの分布）からたまたま選んだ300個ほどのデータから求めた”標準偏差” です。これを $\bar{\sigma}_1$ としましょう。他の1年のデータから求めた”標準偏差” もたまたま選んだ300個ほどデータから求めたもので、これを $\bar{\sigma}_2$ とします。当然 $\bar{\sigma}_{1}\neq \bar{\sigma}_{2}$ となっていることでしょう。どの1年のデータを持ってきても、そこから求められる “標準偏差” はデータセット毎に異なる値になります。
こうして幾つかのバラバラの値の”標準偏差” $\bar{\sigma}_{i}$ が得られました。

まとめると

本当の標準偏差 $\sigma$ は絶対に知りえない
観測から得られた $\bar{\sigma}_{i}$ はばらついている

です。そこで

ばらついている $\bar{\sigma}_{i}$ の平均が本当の標準偏差 $\sigma$ になるようにしよう

という指針を採用します。
すると、STDEVP で求めた標準偏差はこの指針に従わないことが示され、STDEV の定義式で求められる標準偏差こそが新の標準偏差の推定値としてふさわしい（この指針を満たしている）ことが分かります！（細かい計算は適当な教科書を見てください）

結論：

データそのもののバラツキ具合：STDEVP
データの抽出元の真の標準偏差の推定値：STDEV

ということです。

サンプルエクセルシート

Stdevp_Stdev.xls
ダウンロード

真の分布は平均 3（セル B1）、標準偏差 2（セル B2）の正規分布です。
そして、この分布からたまたま選んだ10個（C列からK列）のデータセットが全部で10,000セット（5行目から10004行目）あります。そして各データセットの STDEVP（L列）と STDEV（M列）が計算されています。
STDEVP の平均（セル C2）と STDEV の平均（セル D2）を比較してみてください。STDEV の方が真の分布の推定値にふさわしいことが分かります（もっともっとデータセットを用意すれば STDEV の平均は真の標準偏差に近づいていきます）。

特異値分解入門(3)

May 22, 2010

不正な相関行列 ${\boldsymbol C}$ が与えられると何が起きるが考えてみましょう。
前回の解説の “事実3″ で悲劇が起きます。つまり、

${\small \begin{pmatrix}\lambda_1&0&\cdots&0\\0&\lambda_2&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\lambda_{N}\end{pmatrix}=\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}\begin{pmatrix}\sqrt{\lambda_1}&0&\cdots&0\\0&\sqrt{\lambda_2}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{\lambda_{N}}\end{pmatrix}}$

とする時になって、ハッと気が付くのです…「固有値が負だ！」と。
そしてコンピューターは無情にも「負の平方根は計算できません」というメッセージを残して息絶えます。
ちなみに第1回に載せた不正な相関行列、

${\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}$

の固有値は $-1,2,2$ となります。

この困難を打破してとにかく答えを出したい！そこでいよいよ特異値分解の登場です。

先ずは特異値分解で任意の行列（なんと正方行列でなくてもいい！） ${\boldsymbol A}$ はどのように分解されるか見てみましょう。

${\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}$

となります。とりあえずここでは ${\boldsymbol A}$ として実対称正方行列であるところの相関行列 ${\boldsymbol C}$ だけを対象にして順に説明していきましょう。

行列 ${\boldsymbol C}$ から ${\boldsymbol M}={\boldsymbol C}{\boldsymbol C}^\text{t}$ なる行列を計算します。この行列は自動的に、

正方行列
対称行列
半正定値行列

となるんです！（もちろんもとの行列が実数の行列ならば ${\boldsymbol M}$ も実数）。最後の特徴に注目してください。簡単に言うと、元の行列（ ${\boldsymbol C}$ ）2乗することで負だった固有値があったとしてもそれが正になるということです。
というわけで行列 ${\boldsymbol M}$ は互いに直交する固有ベクトルを持ち、固有値は全て0以上となります。この固有ベクトルを並べた行列が ${\boldsymbol U}$ です。
次に、 ${\boldsymbol N}={\boldsymbol C}^\text{t}{\boldsymbol C}$ なる行列を作ります。この行列は ${\boldsymbol M}$ の場合と全く同様に

正方行列
対称行列
半正定値行列

となります。 ${\boldsymbol N}$ の固有ベクトルを並べた行列が ${\boldsymbol V}$ です。
最後に ${\boldsymbol W}$ は行列 ${\boldsymbol C}$ の固有値の絶対値を対角に並べた対角行列になります。

ここまでの手順で分かるように、各行列 ${\boldsymbol U}$ 、 ${\boldsymbol V}$ 、 ${\boldsymbol W}$ は元の行列 ${\boldsymbol C}$ がどんなものであれ安定的に得ることができます（計算は止まらない！）。前にも述べましたがこれが最も重要な点です。

とにもかくにも行列 ${\boldsymbol C}$ は、

${\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}$

という形に分解されました。次に

${\small {\boldsymbol W}=\begin{pmatrix}|\lambda_1|&0&\cdots&0\\0&|\lambda_2|&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&|\lambda_{N}|\end{pmatrix}=\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}\begin{pmatrix}\sqrt{|\lambda_1|}&0&\cdots&0\\0&\sqrt{|\lambda_2|}&\cdots&0\\\vdots&\vdots&\ddots&\vdots\\0&0&\cdots&\sqrt{|\lambda_{N}|}\end{pmatrix}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime$

とします（対角成分は固有値の絶対値なので安心して平方根がとれる）。すると、

${\boldsymbol C}={\boldsymbol U}{\boldsymbol W}{\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol V}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime({\boldsymbol V}{\boldsymbol W}^\prime)^\text{t}$

となります。

あと一息です

さてここで行列 ${\boldsymbol C}$ が正しい相関行列、つまり半正定値だとどうなるか。その場合は ${\boldsymbol U}={\boldsymbol V}$ となって特異値分解は前回の解説の結果と同じになります。つまり ${\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime$ とすることで、相関行列は ${\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$ となって目標達成です。
しかしそうでない場合はどうでしょう。その場合もやはり同じく ${\boldsymbol Q}={\boldsymbol V}{\boldsymbol W}^\prime$ とすることで手をうちませんか？もちろんこの場合は ${\boldsymbol C}\neq {\boldsymbol Q}{\boldsymbol Q}^\text{t}$ ではありますが（ ${\boldsymbol U}\neq{\boldsymbol V}$ だから）…でも、

安定的に答えが出る
正しい相関行列の場合は正しい答えが出る

という望ましい性質を持ってます。
因みに正しい相関行列ではない場合にも、 ${\boldsymbol Q}{\boldsymbol Q}^\text{t}$ の各成分の符号は元の行列 ${\boldsymbol C}$ のそれに等しくなります（正の相関が負になることはない）。

再び

${\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}$

を例に計算してみましょう。特異値分解の結果は、

${\boldsymbol C}=\underbrace{\begin{pmatrix}0.577&-0.739&-0.348\\-0.577&-0.671&0.466\\0.577&0.068&0.814\end{pmatrix}}_{{\boldsymbol U}}\underbrace{\begin{pmatrix}1&0&0\\0&2&0\\0&0&2\end{pmatrix}}_{{\boldsymbol W}^\prime}\underbrace{\begin{pmatrix}-0.577&0.577&-0.577\\-0.739&-0.671&0.068\\-0.348&0.466&0.814\end{pmatrix}}_{{\boldsymbol V}^\text{t}}$

となります（ ${\boldsymbol U}\neq {\boldsymbol V}$ であることが分かりますか？第1固有ベクトルの符号が逆転しています！こうすることで、固有値の方を正に揃えているという訳です）。そして、

${\boldsymbol Q}=\underbrace{\begin{pmatrix}-0.577&-0.739&-0.348\\0.577&-0.671&0.466\\-0.577&0.068&0.814\end{pmatrix}}_{\boldsymbol V}\underbrace{\begin{pmatrix}1&0&0\\0&\sqrt{2}&0\\0&0&\sqrt{2}\end{pmatrix}}_{{\boldsymbol W}^\prime}}$

となります。ここから ${\boldsymbol Q}{\boldsymbol Q}^\text{t}$ を計算すると、

${\boldsymbol Q}{\boldsymbol Q}^\text{t}=\begin{pmatrix}5/3&1/3&-1/3\\1/3&5/3&1/3\\-1/3&1/3&5/3\end{pmatrix}$

となります（うーん、確かに ${\boldsymbol C}$ とはちょっと違うなぁ。残念）。

以上で相関行列の分解に特異値分解を用いる理由とその方法の解説を終了します。

特異値分解入門(2)

May 21, 2010

ここで相関行列 ${\boldsymbol C}$ が満たしているはずの性質を今一度列挙してみましょう。

成分 $\rho_{ij}$ は実数
対称行列 $\rho_{ij}=\rho_{ji}$
半正定値
対角成分 $\rho_{ii}=1$
対角成分以外 $-1\leq\rho_{ij}\leq 1$

この相関行列 ${\boldsymbol C}$ を、

${\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$

という形に分解するのが目的です（覚えてましたか？）
前回その方法としてコレスキー分解を紹介しましたが、どうも安心して使えなさそうだとわかりました。
次に紹介するのは固有値分解です。これは半正定値実対称行列を分解してくれます。
ではここからは線形代数の勉強です（証明などの詳細は省きます。教科書に必ず載っていますので調べてみてください）。

事実1

任意の実正方行列 ${\boldsymbol A}$ は、ある正方行列 ${\boldsymbol U}$ とこれまた正方行列、ただし対角成分のみしか値がない対角行列 ${\boldsymbol W}$ を使って、

${\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}$

という形に分解される（これを行列のスペクトル分解といいます）。
ここで行列 ${\boldsymbol U}$ は行列 ${\boldsymbol A}$ の固有ベクトルを横に並べたもので、対角行列 ${\boldsymbol W}$ の成分は行列 ${\boldsymbol A}$ の固有値です。

事実2

更に行列 ${\boldsymbol A}$ が対称行列であった場合、行列の固有ベクトルが全て直交する。すると、

${\boldsymbol U}^{-1}={\boldsymbol U}^\text{t}$

となる。というわけで

${\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}$

事実3

更に更に行列 ${\boldsymbol A}$ が半正定値だとすれば、行列の固有値は全て0以上。つまり行列 ${\boldsymbol W}$ は、

と出来るではないか！すると、

${\boldsymbol A}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^{-1}={\boldsymbol U}{\boldsymbol W}{\boldsymbol U}^\text{t}={\boldsymbol U}{\boldsymbol W}^\prime {\boldsymbol W}^\prime {\boldsymbol U}^\text{t}$

となります。ここで $({\boldsymbol A}{\boldsymbol B})^\text{t}={\boldsymbol B}^\text{t}{\boldsymbol A}^\text{t}$ という事実（積の順番が入れ替わる）と、対角行列は転置しても変化しないこと（この場合 $({\boldsymbol W}^\prime)^\text{t}={\boldsymbol W}^\prime$ ）に注意すると、

やった！！ ${\boldsymbol U}{\boldsymbol W}^\prime = {\boldsymbol Q}$ とおけば目標達成じゃないか！

${\boldsymbol A}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$

しかし…
実務上（ここがポイント）、相関行列がいつもいつもタチの良い形であるとは限らないのです！もちろん実対称行列であることはどんな場合も間違いないですが（複素数の相関や、 $100\times 234$ の相関行列を作れるものなら作ってみろ）、 半正定値を満たさない場合は往々にしてありえるのです。例えば…

データに欠損があって、欠損データを適当に埋める
乱数列のうち、いくつかのペアについては相関が（観測とは無関係に）与えられている
測定誤差

などの場合、相関係数に矛盾が出ることがあるのです。相関係数が矛盾するとはどういうことか、ひとつ最も大げさな例を示してみましょう。
以下の3変数の相関行列を見てください。

${\boldsymbol C}=\begin{pmatrix}1&1&-1\\1&1&1\\-1&1&1\end{pmatrix}$

実数であるとか、対称であるとか相関行列の性質は確かに満たしているんですが…何がおかしいか分かりますか？

第1変数と第2変数の相関係数 $\rho_{12}$ は 1。第1変数と第2変数は完全相関で全く同じ値になる
第2変数と第3変数の相関係数 $\rho_{23}$ は 1。第2変数と第3変数は完全相関で全く同じ値になる

ってことは必然的に第1変数と第3変数は同じ値になるはずです。しかし与えられた相関係数 $\rho_{13}$ は -1 になっている！つまり矛盾しているわけです。
要するに相関係数がありえない組み合わせになっている場合がありえるのです（さすがに例のような大げさなのはないでしょうけど）。このとき、相関行列は半正定値ではなくなってコレスキー分解が機能しなくなります。

実務上、相関行列が与えられたら（それに矛盾があろうとなかろうと）計算がストップしてしまうことは望ましくありません。とにかく安定して何らかの答えを出すことが求められます。もちろん正しい相関行列（半正定値実対称行列）は正しく分解されることが前提です。

そこで NtRand では相関行列の分解に特異値分解(Singular Value Decomposition，SVD)という手法を採用しています。

特異値分解入門(1)

May 20, 2010

互いに相関を持った正規乱数のセット（多変量正規乱数）を生成すること、これが金融工学のモンテカルロシミュレーションにおいて最も重要な技術となります。 $N$ 系列の正規乱数列を生成する場合を考えます。これらの乱数列の相関構造を記述するのが相関行列と呼ばれる $N\times N$ の正方行列 ${\boldsymbol C}$ です。最終的に

${\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$

という形に分解することが目標となります（ここで t は転置行列です）。分解して得られた行列 ${\boldsymbol Q}$ を使って実際に相関乱数を生成する方法については、これも一般的な数値計算の教科書に載っていますのでそちらをご覧ください。

すぐ後で説明するように、この分解にはコレスキー分解を使えと、多くの教科書で教えていますが、実務上コレスキー分解ではうまくいかない場合があります。
NtRand あるいは Numerical Technologies 社の製品では相関行列の分解に特異値分解という方法を採用しています。
これから3回にわたって、なぜコレスキー分解ではダメなのか、そしてなぜ特異値分解だとうまくいくのかを説明します。

相関行列について

相関行列 ${\boldsymbol C}$ は、その ij 成分が確率変数 $X_{i}$ と $X_{j}$ との相関係数 $\rho_{ij}$ である行列です。この定義から対角成分 $\rho_{ii}$ は確率変数 $X_{i}$ と自分自身との相関、つまり 1 となります。また、 $\rho_{ij}=\rho_{ji}$ も定義から明らかです。
つまりこの行列は、

$N\times N$ 正方行列
成分は実数
対角成分は 1
対角成分以外は $-1\leq\rho_{ij}\leq 1$
対称（ $\rho_{ij}=\rho_{ji}$ ）

という性質を持っています。
更に相関係数の定義から、この行列は半正定値であることが示されます（半正定値って何だ？それは行列の固有値が全て0以上であること。固有値って何かって？それは勉強してください）。つまり相関行列はその定義から半正定値実対称行列であるということです。

コレスキー分解 (Cholesky decomposition)

${\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$

という分解を行う技として、どの教科書にも書いてあるのがコレスキー分解という方法です。どの教科書にも載っているので、ここでは説明しません。
この手法は正定値実対称行列を目的の形に分解してくれるので、正定値実対称行列である相関行列に対して使えば目標達成！メデタシ、メデタシとなる筋書きです。
しかし…
実はコレスキー分解は容易に失敗します。その代表的な状況が観測データが少ない場合です。確率変数の数が100個で、その観測データの数が100個未満の場合は、定義に従った正しい相関行列であってもコレスキー分解は失敗します。なぜなら、この状況で作られる相関行列は必然的に半正定値行列になるからです。コレスキー分解は正定値行列専用です。
それでは、少なくとも正しく作られた相関行列では間違いなく

${\boldsymbol C}={\boldsymbol Q}{\boldsymbol Q}^\text{t}$

と分解してくれるより強力な方法はないでしょうか？

火星に猫はいるか？

男の子？女の子？

メルセンヌ素数

サンクトペテルブルグのパラドックス