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

May 25, 2010

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

「サイコロを転がして出た目が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} は、
{\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}}={\boldsymbol W}^\prime {\boldsymbol W}^\prime
と出来るではないか！すると、
{\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 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 U}{\boldsymbol W}^\prime ({\boldsymbol U}{\boldsymbol W}^\prime)^\text{t}
やった！！ {\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}
と分解してくれるより強力な方法はないでしょうか？
 

コンカツの数学

May 19, 2010

モテモテで我が世の春を謳歌していたあなたは、いよいよ真剣に結婚を考えることにしました。しかし焦って貧乏くじを引きたくない！できるだけいい条件の人と結婚したい！
あなたは勝負はこの1年と決めて、ある結婚相談所に登録することにしました。そこはかなり値段が張るだけはあって「今後1年で100人のお相手の紹介を保証」とのこと。ただし…
  • (1) 前もって相手のプロフィールなどは見せられない（条件などは会うまで分からない）
  • (2) 1度に1人しか会えない
  • (3) 交際を断った場合に次の人を紹介する
  • (4) 一度断った相手とは二度と会えない（キープは不可）
  • (5) こちらが OK すれば相手は必ず応じる。そしてゴールイン！
という条件付きです。
この人だ！と思ってももしかしたら次にもっといい人が現れるかも…でもそうこうしているうちに1番の人を見逃してしまっているかも…これは綿密なる戦略が必要だ！

そこであなたは先ず勇気を出して最初の何人かを見逃して（観察して）、その後「これまでに見逃した人より高条件の人」が現れたら問答無用でその人をゲットする、という戦略を採用することにしました。

この戦略だと、見逃す人数が少なすぎると、高条件の人を見逃す確率は減るけども低条件の人を選んでしまう確率も増える。見逃す人が多すぎると、最高条件の人を見逃してしまうかもしれない…では何人を見逃すことで最高の相手をゲットする確率が一番高くなるでしょうか？

結論を先に言うと、最初の37人を見逃すのが最良の戦略で、その場合最高の相手をゲットする確率は約37.1%となります。

ここから先は数学の話

N 人の見合い相手（今の例では N=100）のうち、見逃す人数を s\;(0\leq s&lt N) とする。
本当の No.1 の相手が j\;(j=1,2,\cdots,N) 番目にいる場合で、しかもその人を選ぶ確率を P_{j} と書くとしよう。
P_{j}を考えてみよう。これは

（j 番目に最高の相手がいる確率）×（そこに達する確率）

と定義される。
前者の確率は 1/N となる（完全にデタラメなので一様分布）。
後者の確率は、運命の j 番目の相手に至るまでに出会う（j-1）人のうちの最高の人が、最初の s 人に含まれていればいいというのは分かるだろうか？ もし（j-1）人のうちの最高の人が最初の s 人に含まれていなければその人に出会った時点で婚活は終了し、運命の人には出会えないことになるからである。 逆にこの人が最初の s 人に含まれていれば、観察期間が終了して運命の j 番目の相手に出会うまでは見逃した人より高条件の人には出会わないことになる。 というわけで、後者の確率は s/(j-1) である。結局 P_{j} は、
P_{j}=\frac{1}{N}\times\frac{s}{j-1}\;(j=s+1,\;s+2,\;\cdots\;N)
となる（j が 1 から s までの場合は最高の人を見逃してしまった！ということで出逢う確率は 0）。
最高の相手に出会う確率は P_{j} の和、つまり
P=\sum_{j=s+1}^{N}P_{j}=\frac{s}{N}\sum_{j=s+1}^{N}\frac{1}{j-1}
となる。あとはこれを最大にする s を見つければよいことになる！

ちなみにこの 37% という何ともキレの悪い数字は一体何かというと、
\frac{1}{\text{e}}=0.367879\cdots
です。何とこんなところに自然対数の底 \text{e}=2.71828\cdots が！！
 

黄金のフラフープ

May 18, 2010

この一様乱数（連続）を利用したゲームは
http://beetama.blog14.fc2.com/blog-date-200807.html
に掲載されていたものを参考にしました。ペコリ。

ルール

1箇所に赤い丸印（基準点）がつけてあるリングがあります。リングの任意の2箇所にボタンを使って切り込みを入れます。その結果リングは2つに分割されますが、そのうち赤い丸印がない方があなたの取り分、もう一方が闇のオーナーの取り分となります。

  • あなたの取り分が大きい場合：あなたの勝利。闇のオーナーとの取り分の差額が手に入る。
  • あなたの取り分が小さい場合：あなたの負け。闇のオーナーとの取り分の差額を支払う。
要するに基準点、切り込み1、切り込み2の3点を円周上にランダムに配置しているだけなのです。その位置は一様分布（連続）に従っていて、範囲は 0\sim 360[度]になっているというわけ。 3点をデタラメにばら撒いているだけなので、どっちに基準点が入るかは5分5分のような気がする。ならば勝敗は5分5分なんでしょうか？

[スタート]ボタンをクリックするとリングが回転します。その後、同じボタン（[最初の切り込み]となっている）をクリックすると赤い三角形の箇所に切り込みが入り、更に同じボタン（[次の切り込み]となっている）をクリックすると2箇所目に切り込みが入ってリングは停止します。 あなたの取り分が緑色で示されます。あなたは見事勝利することができるでしょうか？勝率や成績を調べてみてくださいね（果たして5分5分？）。