元バイオ系

元バイオウェット系がデータサイエンスやらを勉強していくブログ。 基本自分用のまとめ。

2次形式の最大値と最小値

数学 Python

前回の記事で、固有値と固有ベクトルについて書いた。

hotoke-x.hatenablog.com

この記事中にある正方行列の対角化を式変形すると、行列 \displaystyle Aのスペクトル分解（固有値分解）が得られ、これを使って2次形式の標準形が得られる。

正方行列のスペクトル分解（固有値分解）

正方行列 \displaystyle Aの対角化は以下のような式だった。

$$ \begin{align} \boldsymbol{P^{-1}AP} &= \boldsymbol{D} \\ &= \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{align} $$

\displaystyle Pは固有ベクトル、 \displaystyle Dは対応する固有値を対角成分に並べたものである。
これを式変形して、 \displaystyle Aのスペクトル分解を得る。

$$ \begin{align} \boldsymbol{A} &= \boldsymbol{P D P^{-1}} \end{align} $$

2次形式の標準形

正方行列のスペクトル分解を使うと、2次形式の標準形が得られる（2次形式と2次関数は別物なので注意）*。

まず、2次形式

$$ \begin{align} \boldsymbol{x^{T} A x} \end{align} $$

を考える。次に、ベクトル \displaystyle \boldsymbol{x}を固有ベクトルを基底とする空間 \displaystyle \boldsymbol{P}に投げる。

$$ \begin{align} \boldsymbol{t} &= \boldsymbol{P x} \\ \boldsymbol{x} &= \boldsymbol{P^{T} t} \end{align} $$

これを2次形式の式に代入すれば

$$ \begin{align} \boldsymbol{x^{T} A x} &= \boldsymbol{\left(P^{T} t \right)^{T} A P^{T} t} \\ &= \boldsymbol{t^{T} P A P^{T} t} \\ &= \boldsymbol{t^{T} D t} \end{align} $$

となる。ただし

$$ \begin{align} \boldsymbol{P^{-1}} &= \boldsymbol{P^{T}} \\ \boldsymbol{P^{-1}AP} &=\boldsymbol{PAP^{-1}} = \boldsymbol{D} \end{align} $$

を利用した。これを2次形式の標準形と呼ぶ。ただ、固有値と固有ベクトルを使って2次形式を書き直しただけである。

二次形式の最大値、最小値

\displaystyle \boldsymbol{A}が半正定値対称行列なら、2次形式の最大値、最小値はそれぞれ \displaystyle \boldsymbol{A}の最大固有値、最小固有値に対応する。標準形を得ると、このように話がシンプルになる。

少々雑だったが、ここまでがよくある書籍の解説である。これですっきり理解できるならそれでよいが、なんで最大固有値、最小固有値が最大値、最小値に対応するのかイメージできないかもしれない。

2次形式の標準形の幾何学的解釈

これなら分かる最適化数学（金谷健一）、例題1.25、1.27より \displaystyle x^{2}+y^{2}=1の時

$$ \begin{align} f=6x^{2}+4xy+3y^{2} \end{align} $$

の最大値、最小値を考える。これを標準形に直すと

$$ \begin{align} f=2x' ^2 + 7y' ^2 \end{align} $$ となる。 \displaystyle x'、 \displaystyle y'は固有ベクトル空間 \displaystyle \boldsymbol{P}で変換後の \displaystyle x、 \displaystyle yで

$$ \begin{align} \boldsymbol{t} &= \boldsymbol{P x} \\ \boldsymbol{x} &= \left(\begin{array}{c} x \\ y \end{array}\right), \boldsymbol{t} = \left(\begin{array}{c} x' \\ y' \end{array}\right) \\ \end{align} $$

である。係数2と7は固有値になっている。

また、 \displaystyle x^{2}+y^{2}=1と x'^{2}+y'^{2}=1は同値である。

改めて標準形の式を眺めてみると、幾何学的には放物面に対応していることがわかる。実際にプロットしてみると

f:id:hotoke-X:20190103202346p:plain

となり、等高線が楕円になっていることがわかる。 \displaystyle \boldsymbol{x'^{2}+y'^{2}}=1に制限すれば、 \displaystyle x'=1のとき最小値2（最小固有値）、 \displaystyle y'=1のとき最大値7（最大固有値）になることがわかる。

*2次関数は2次以下の項からなる関数。2次形式は2次の項のみからなる関数。

参考書籍

  • 金谷健一（2005）「これなら分かる最適化数学」共立出版株式会社

固有値と固有ベクトル

数学

固有値と固有ベクトルについてメモ。
何をしているかはわかっているけど、何が起きているかがわかっていなかった。

  • 「固有」ってなんだよと思っている
  • 行列の掛け算で何が起きてるのかイメージがわかない

そんな人たちは参考になるかも。

数弱なので間違いがあれば教えてくれると助かります。

固有値と固有ベクトルについておさらい

行列 {\displaystyle A} に対して $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$

が成り立つ{\displaystyle \boldsymbol{0}}でないベクトル{\displaystyle \boldsymbol{u}}を行列{\displaystyle A}の固有ベクトル、{\displaystyle \lambda}をその固有値という。
{\displaystyle A}が{\displaystyle n \times n}対称行列なら、{\displaystyle A}は{\displaystyle n}個の固有値{\displaystyle \lambda_1,\ldots,\lambda_n}を持つ。また、{\displaystyle \boldsymbol{u}}の各成分間は互いに直交する。

ここで、初めの式をちょっと変形する $$ \begin{equation} \left(\lambda \boldsymbol{I} - \boldsymbol{A} \right) \boldsymbol{u} = \boldsymbol{0} \end{equation} $$ 定義より{\displaystyle \boldsymbol{u}} \neq \boldsymbol{0}より $$ \begin{equation} | \lambda \boldsymbol{I} - \boldsymbol{A}| = \boldsymbol{0} \end{equation} $$

この式を固有方程式（特性方程式）と呼ぶ。

つまりなんなのか

行列{\displaystyle A}を左からかけることを、行列{\displaystyle A}を作用させるなんて物理では言うらしい。
つまり、ベクトルを伸ばしたり、縮めたり、回転したりして別のベクトルへ変換することを意味している。

ここで改めて式を眺めてみる。 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ この式が意味するところは、「行列{\displaystyle A}による{\displaystyle \boldsymbol{u}}の変換は、結局{\displaystyle \boldsymbol{u}}を{\displaystyle \lambda}倍することと同じ」ということ。
別の言い方をすれば、「行列による変換をしても、長さが{\displaystyle \lambda}倍変化するだけで向きが変わらないベクトル{\displaystyle \boldsymbol{u}}がある」ということで、そんなベクトルを固有ベクトルと呼んでいる。

基底変換

今、デカルト座標系にベクトル $$ \begin{align} \boldsymbol{v} = \left( \begin{array}{c} x \\ y \end{array} \right) \end{align} $$ が存在して、座標変換によって新たな表現を得ることを考える。
これはデカルト座標での表現なので、基底ベクトル{\displaystyle \boldsymbol{e}_x}, {\displaystyle \boldsymbol{e}_y}を使って大げさに書けば $$ \begin{align} \boldsymbol{v} &= x\boldsymbol{e}_x + y\boldsymbol{e}_y \\ &= x\left( \begin{array}{c} 1 \\ 0 \end{array} \right) +y\left( \begin{array}{c} 0 \\ 1 \end{array} \right) \\ &= \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right) \left(\begin{array}{c} x \\ y \end{array}\right) \end{align} $$ となる。
この基底ベクトルを $$ \begin{align} \left(\begin{array}{cc} 1 & 0 \\ 0 & 1 \end{array}\right)\left(\begin{array}{c} x \\ y \end{array}\right)=\left(\begin{array}{cc} p_x & q_x \\ p_y & q_y \end{array}\right)\left(\begin{array}{c} x' \\ y' \end{array}\right) \end{align} $$ となるような別の基底ベクトル $$ \begin{align} \boldsymbol{p}=\left(\begin{array}{c} p_x \\ p_y \end{array}\right), \;\; \end{align} \boldsymbol{q}=\left(\begin{array}{c} q_x \\ q_y \end{array}\right) $$ で置き換えて、xy座標（デカルト座標）以外の任意のpq座標で表現できることがわかる。
ここでまた、最初の式をまた眺めてみる。 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ この式の正方行列{\displaystyle A}は、ベクトルの変換意味しているのだった（ベクトル空間の線形写像）。
説明のために、 $$ \begin{equation} \boldsymbol{A} = \left(\begin{array}{cc} a_{11} & a_{12} \\ a_{21} & a_{22} \end{array}\right), \;\; \boldsymbol{u} = \left(\begin{array}{c} u_1 \\ u_2 \end{array}\right) \end{equation} $$ とすると、 $$ \begin{align} \boldsymbol{A u}=\left(\begin{array}{cc} a_{11}\\ a_{21} \end{array}\right)u_1 +\left(\begin{array}{cc} a_{21} \\ a_{22} \end{array}\right)u_2 \end{align} $$ と書けることから、行列{\displaystyle A}の列成分は変換後の空間での基底ベクトルを表していたことがわかる。 つまり、{\displaystyle (u_1, u_2)}という数字をそのまま行列{\displaystyle A}による変換後の空間に渡すとどうなるか計算していただけのことだった。

正方行列の対角化

さて、行列{\displaystyle A}による変換前後で固有ベクトル{\displaystyle \boldsymbol{u}}はスケールが固有値{\displaystyle \lambda}倍変わるだけで向きが変わらないのであった。
ということは、

  1. 行列{\displaystyle A}による固有ベクトルのスケール変化だけを事前に調べておく（行列Aの固有ベクトルと固有値を調べておく）
  2. 変換後の空間を、固有ベクトルを基底として表現する

という手順で行列{\displaystyle \boldsymbol{A}}による変換は簡単な式に置き換えられそうだ。

{\displaystyle n\times n}正方行列{\displaystyle A}によってn次元ベクトル{\displaystyle \boldsymbol{s}}が別のn次元ベクトル{\displaystyle \boldsymbol{t}}に変換される、 $$ \boldsymbol{t} = \boldsymbol{A s} $$ という変換を考える。ここに、行列{\displaystyle A}の右側に固有ベクトルを並べた行列 $$ \begin{align} P&=\left(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_n\right) \\ \boldsymbol{p}_1 &= \left(\begin{array}{c} p_{11} \\ \vdots \\ p_{1n} \end{array}\right), \;\; \boldsymbol{p}_2 = \left(\begin{array}{c} p_{21} \\ \vdots \\ p_{2n} \end{array}\right), \;\; \cdots, \;\; \boldsymbol{p}_n = \left(\begin{array}{c} p_{n1} \\ \vdots \\ p_{nn} \end{array}\right) \end{align} $$ をねじ込んでみると $$ \boldsymbol{t} = \boldsymbol{(AP)(P^{-1}s)} $$ となる。{\displaystyle \boldsymbol{AP}}の部分は、一番最初に登場した式 $$ \begin{equation} \boldsymbol{A u} = \lambda \boldsymbol{u} \end{equation} $$ の左辺に対応している。異なるのは、固有ベクトルをひとつずつ扱っているのではなく、行列として並べてしまったことくらいである。じゃあ右辺はどうなってるのと思うかもしれないが、固有値と対応する固有ベクトルの積が以下のように並んでいるだけだ。 $$ \begin{equation} \boldsymbol{A P} = A \left(\boldsymbol{p}_1, \ldots, \boldsymbol{p}_n\right) = \left(\lambda_1 \boldsymbol{p}_1, \ldots, \lambda_n\boldsymbol{p}_n\right) \end{equation} $$

話を戻す。
$$ \boldsymbol{t} = \boldsymbol{(AP)(P^{-1}s)} $$ の式をよく見てみると、{\displaystyle \boldsymbol{s}\longrightarrow\boldsymbol{t}}とする変換だったのに、{\displaystyle \boldsymbol{P^{-1}s}\longrightarrow\boldsymbol{t}}のように見えてしまう。{\displaystyle \boldsymbol{AP}}を使いたいために格好が悪くなってしまった。

いっそのこと、{\displaystyle \boldsymbol{P^{-1}s}\longrightarrow\boldsymbol{P^{-1}t}}と考えてしまえばよさそうだ。つまり、{\displaystyle \boldsymbol{s}} と{\displaystyle \boldsymbol{t}} をそれぞれ{\displaystyle \boldsymbol{P^{-1}}} が基底になっている空間に投げてしまってから考えるということだ。もし、元の座標系での値が知りたかったら{\displaystyle \boldsymbol{P}}を掛ければ元通りなので問題はない。 式で書くと

$$ \begin{equation} \boldsymbol{P^{-1} t} = \boldsymbol{(P^{-1}AP)(P^{-1}s)} \end{equation} $$

となる（有名な式らしい...知らんけど）。

なんだか大事になっているように見えるが、実はそうではない。\displaystyle \boldsymbol{P^{-1}AP}の部分が実は、固有値を対角成分に並べた対角行列になるからだ（\displaystyle \boldsymbol{P^{-1}AP}を対角行列にする操作を対角化と呼ぶ）。

$$ \begin{align} \boldsymbol{(P^{-1}AP)} &= \boldsymbol{D} \\ &= \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{align} $$

つまり、行列{\displaystyle A}によるベクトルの変換は、{\displaystyle \boldsymbol{P^{-1}}}が基底となる空間へ投げてしまえば、そこでは基底ベクトルの方向へ固有値倍しただけだとみなせるということである。

「固有ベクトルを基底とした空間へ投げるんじゃないの？」

と思われるかもしれないがこれは書き方の問題で、結局

$$ \begin{align} \boldsymbol{P^{-1}AP} =\boldsymbol{PAP^{-1}} = \boldsymbol{D} \\ \end{align} $$

なので、

$$ \begin{equation} \boldsymbol{P^{-1} t} = \boldsymbol{(P^{-1}AP)(P^{-1}s)} \end{equation} $$

でも

$$ \begin{equation} \boldsymbol{Pt} = \boldsymbol{(PAP^{-1})(Ps)} \end{equation} $$

でも同じことである。「固有ベクトルを基底とした空間に投げてから考える」という意味では後者の方がスッキリしているのかもしれない。

【数理統計】Fisher情報量とクラメール・ラオ下限（CR下限）

統計

Fisher情報量って何なのか理解したかっただけ。
Fisher情報量がわかれば、クラメール・ラオ下限（不偏推定量の分散の下限）がわかる。
備忘録として簡単にまとめておく。

前提として、不偏推定量についてだけ考えるものとする。
不偏推定量に関してはこちら

hotoke-x.hatenablog.com

推定量が不偏性を持てば、その推定量の最適性・適切性の判定が容易になる。 互いに独立な標本{\displaystyle　X_1, \ldots, X_n}の同時密度関数が{\displaystyle f(x|\theta)}で与えられるとする。すなわち $$ \begin{align} X_i \in \boldsymbol{X} &, \quad i \in \left\{1, \ldots, n \right\} \\ f(\boldsymbol{x}|\theta) := \prod_{i=1}^n p(x_i; \theta) &,\quad \boldsymbol{x} = \left\{x_1, \ldots, x_n \right\} \end{align} $$

このとき、 $$ \begin{align} \forall \theta \in \Theta,\qquad \boldsymbol{E} \left(\hat \theta_n \right) \end{align} $$

を満足する。これを{\displaystyle \theta}について微分 $$ \begin{equation} \frac{\partial}{\partial \theta} \boldsymbol{E} \left(\hat \theta_n \right) = 1 \end{equation} $$

これを展開 $$ \begin{equation} \frac{\partial}{\partial \theta} \int f(\boldsymbol{x}|\theta) ~ \hat \theta_n ~ \mathrm{d}\boldsymbol{x} = \int \hat \theta_n \frac{\partial}{\partial \theta} f(\boldsymbol{x}|\theta) \mathrm{d}\boldsymbol{x} \label{eq:fisher1} \end{equation} $$

ここで、 $$ \begin{align} \frac{\partial}{\partial \theta} \log{L} &= \frac{1}{L} \frac{\partial L }{\partial \theta} \\ \frac{\partial L}{\partial \theta} &= L \frac{\partial}{\partial \theta} \log{L} \end{align} $$

より

$$ \begin{align} \mathrm{R.H.S} &= \int \hat \theta_n \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} f(\boldsymbol{x}|\theta) \mathrm{d} x \\ &= \boldsymbol{E} \left[\hat \theta_n \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right] \\ &= Covariance \left(\hat \theta_n, \frac{\partial \log{f(\boldsymbol{x} | \theta)} }{\partial \theta} \right) \\ &\leq \sqrt{V \left(\hat \theta_n\right)} \sqrt{I_n\left(\theta \right)} \label{eq:fisherinfo} \end{align} $$

ここで、{\displaystyle I_n\left(\theta \right)}はFisher情報量と呼ばれる量で、クラメール・ラオの下限（CR下限）と密接な関係がある（後述）。なお、

$$ \begin{align} Covariance \left(\hat \theta_n, \frac{\partial \log{f(\boldsymbol{x} | \theta)} }{\partial \theta} \right) \leq \sqrt{V \left(\hat \theta_n\right)} \sqrt{I_n\left(\theta \right)} \end{align} $$ ではコーシー・シュワルツ不等式を利用した。また、

$$ \begin{align} I_n\left(\theta \right) &= V \left(\frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) \\ &= \boldsymbol{E} \left[\left( \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta}\right) ^2 \right] \end{align} $$

特に{\displaystyle n}個の標本が{\displaystyle i.i.d}なら $$ \begin{align} I_n \left( \theta \right) &= \boldsymbol{E} \left[ \left( \frac{\partial}{\partial \theta} \sum_{i=1}^{n} \log{f(\boldsymbol{x}|\theta)} \right) ^2 \right] \\ &= n \boldsymbol{E} \left[\left(\frac{\partial}{\partial \theta} \log{f(\boldsymbol{x_1}|\theta)} \right) ^2 \right] \\ &= nI_1 \left( \theta \right) \end{align} $$

ここで $$ \begin{equation} \frac{\partial}{\partial \theta} \boldsymbol{E} \left(\hat \theta_n \right) = 1 \end{equation} $$

を思い出せば $$ \begin{align} Covariance \left(\hat \theta \left(\boldsymbol{x} \right), \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) &= 1 \\ Covariance \left(\hat \theta \left(\boldsymbol{x} \right), \frac{\partial \log{f(\boldsymbol{x}|\theta)}}{\partial \theta} \right) &\leq \sqrt{V \left(\hat \theta_n \right)} \sqrt{I_n \left(\theta \right)} \end{align} $$

より $$ \begin{align} 1 &\leq \sqrt{I_n \left(\theta \right)} \\ V \left(\hat \theta \right) &\geq \frac{1}{nI_1 \left(\theta \right)} \label{eq:CR} \end{align} $$

となる。以上より、フィッシャー情報量{\displaystyle I_1 (\theta)}によって推定量の分散の下限が決まることがわかる（推定量としての良さがわかる）。この下限をクラメール・ラオ下限（CR下限）と呼ぶ。

【Julia】Julia-1.0.0のProjectを試す【仮想環境？】

Julia

Julia1.0.0が遂にリリースされました。

f:id:hotoke-X:20180810011534p:plain

仮想環境っぽいものを作って環境の切り替えができるようなので、試してみます。

Projectの作成

Anaconda Pythonではconda createコマンドで仮想環境を作れました。

Julia1.0.0ではProjectと呼ぶようです。

Anacondaでは

conda create -n envname

とする必要がありましたが、Juliaではプロジェクトのディレクトリに入って

activate .

とするようです。

その後、

julia>]
(projectname)pkg> add hogehoge

で、プロジェクトローカルな環境にhogehogeパッケージが入ります。

その後バックスペースでpkg環境を抜けて同フォルダで作業すれば良いっぽい。

ただ、Julia1.0.0がリリースされたのが記事作成時点で昨日の今日なので、対応が追い付いているパッケージを探すのが結構大変です（2018/08/10現在）。

私は情弱なのでおとなしく修正を待つことにします。

(ちなみに現状遊ぶなら0.6が無難)

文献検索、管理、メモを1つの画面で完結させる

研究ライフハック

文献管理ソフトの紹介です。

以前にPaperpileが凄いって紹介をしました。 hotoke-x.hatenablog.com

個人的には好きなんですが、wordとの相性が悪い...

もっと良いのはないかと思ってたらReadcubeがありました。

数か月使ったので良い部分を簡単に紹介していきます。

※例示しているのはオープンアクセスな論文です。

Readcubeの凄いところ

  1. 複数の検索エンジンで同時検索
  2. 検索結果が自動保存される（そして簡単に消せる）
  3. 毎年の引用数がわかる
  4. 引用されている文献へ容易に飛べる
  5. どんな風に引用されているのか見れる
  6. Recommendationsが強力

こんな感じで全体的に便利すぎ。
全部1つの画面で完結します。
関連論文が次から次へ見つかって大変（嬉しい）。

使い方

ユーザー登録は済んでいるものとします。
また、ウェブ版での説明です（デスクトップ版もあります）。

基本画面（青い四角は身バレ防止です。）

f:id:hotoke-X:20180623035811p:plain

画面にNatureやらScienceやらGoogle Scholarやらが表示されてますが、こっから同時検索してくれます。 なんかもうこれだけで凄い...。

f:id:hotoke-X:20180623035441p:plain

検索してみた結果。 アブストが表示される。 オープンアクセスだったり、学校などジャーナルを購読している機関からアクセスするとFigureまで表示してくれます。

f:id:hotoke-X:20180623035625p:plain

検索結果の図をクリックすると拡大表示。この段階ではまだ文献をダウンロードしたりしていません。 あくまでも検索結果です。

「Add to Library」をクリックするとインポートしてくれます。
（supplementがあればそれも同時にインポートする。すごい...）

f:id:hotoke-X:20180623035836p:plain

そして、画面左側には検索結果が自動保存されます。便利！ 青い部分には自分でタグ付けした文献が入っています。 一つの文献に複数のタグをもたせることもできます。

ここからが特にお気に入りの機能。

f:id:hotoke-X:20180623040308p:plain ライブラリで文献をクリックすると、情報が表示されます。 アブストがここでざっと確認出来ます。

で、上のグラフのマークをクリックすると引用情報が見れます。 f:id:hotoke-X:20180623040530p:plain

毎年どれくらい引用されているのか、どの論文に引用されているのかを見ることができます。

f:id:hotoke-X:20180623040724p:plain

どんなふうに引用されているかも見ることができる。凄すぎて拍手。

その他

デスクトップ版ではwordへのcitation追加機能も付いてきます。
もちろん文献を読みながらハイライトしたり、ノートをとることも可能です。
デスクトップアプリが重いことが唯一残念。

Recommendationsをクリックすると、ライブラリの文献を元にお勧め論文を検索して表示してくれます。
しかもタグごとのお勧めも表示可能。
関連論文がめっちゃ見つかります。

JuliaでPRML：多項式曲線フィッティング

Julia PRML

Juliaの勉強も兼ねてPRMLの多項式フィッティングを実装してみた。

誤差関数書いたけど、今回は解析的にパラメータが出てくるので結局使っていない。
GAとかMCMCとかでパラメータ推定して性能見てみるのもいいかも。
計画行列をうまく使えるようになりたい。

【数理統計】ラオーブラックウェル化（Rao-Blackwellization）

統計

理解したいのはこっちでした。

本質的には同じ問題だけど。

ラオーブラックウェル化（Rao-Blackwellization）は、条件付き期待値を利用して期待値を計算する方法です。

サンプルが独立なら、ブラックウェルーラオの定理より直接期待値を計算するより精度が良くなる事が保証されます（不思議）。

名前の順番がひっくり返っているのは知りません（書籍でこう書いてあったからそのまま書いてます）

ブラックウェルーラオの定理の証明はこちら

hotoke-x.hatenablog.com

ラオーブラックウェル化（Rao-Blackwellization）

確率変数 \boldsymbol{x}を二つのブロックに分割し、 \boldsymbol{x} = (\boldsymbol{x_1}, \boldsymbol{x_2}) (x_i \in \chi_i)と表す。このとき、 x_1の統計量の期待値

f:id:hotoke-X:20180322193210p:plain

を求めることを考えます。

\pi (\boldsymbol{x_1}, \boldsymbol{x_2})からサンプリングされた \left( \boldsymbol{x^{(1)}}, \ldots, \boldsymbol{x^{(t_0 + T)}} \right)を用いれば、上式の期待値は以下のように計算できます。

\displaystyle \hat I = \frac{1}{T} \sum_{t=1}^{T} g(\boldsymbol{x_1^{(t_0 + T)}})

ここで、 \boldsymbol{x^{(t)}} = \left(\boldsymbol{x_1^{(t)}}, \boldsymbol{x_2^{(t)}} \right)であることに注意。

さらに \boldsymbol{x_2}の条件付き分布で同時分布を分解すると

f:id:hotoke-X:20180322194353p:plain

となるので、この式の真ん中部分

\displaystyle h(\boldsymbol{x_2}) = \boldsymbol{\mathrm{E}} \left[g(\boldsymbol{x_1}) | \boldsymbol{x_2} \right] = \int_{\chi_1} g(\boldsymbol{x_1})\pi(\boldsymbol{x_1} | \boldsymbol{x_2}) \mathrm{d}\boldsymbol{x_1}

を解析的に求められれば、

\displaystyle \hat I_{RB} = \frac{1}{T} \sum_{t=1}^{T} h(\boldsymbol{x_2}^{(t_0 + t)})

によって計算できる。

なんだ面倒くさいと思うが、実はラオーブラックウェルの定理から

\displaystyle \boldsymbol{\mathrm{V}} (\hat I) \geq \boldsymbol{\mathrm{V}} (\hat I_{RB})

が成立し、より良い推定値になり得ることがわかる。

参考書籍

  1. 応用を目指す数理統計学（国友直人）
  2. ベイズ計算統計学（古澄英男）

【数理統計】ブラックウェルーラオ（Blackwell-Rao）の定理

統計

はい。

私がちゃんと理解したかった定理です。

復習したら見た瞬間理解したのですが、一応メモっておきます。

初学者向けの説明がなかなかネットに落ちていなかった（気がする）ので。

初めに行っておくと、この定理は

「十分統計量で条件づけた不偏推定量の分散は、他の不偏推定量の分散より大きくなることはない。」

って定理です（違ったらコメントください）。

十分統計量に依存した不偏推定量を使っとけばとりあえず無難ですよってことですな。

数学アレルギーでも、この程度の認識は持っておいた方が良いでしょう。

十分統計量についてはこちら hotoke-x.hatenablog.com

不偏推定量についてはこちら hotoke-x.hatenablog.com

証明

証明は結構簡単です。

今、未知母数の推定量は不変推定量しか考えないとします。このとき、標本 \boldsymbol{X}について \thetaの不偏推定量 \phi (\boldsymbol{X})、十分統計量 T(\boldsymbol{X})を考えます。この時十分統計量で条件付けした統計量を \varphi (t) = \boldsymbol{\mathrm{E}} \left[\phi (\boldsymbol{X}) | T = t \right]とすると

\boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \theta \right)^2 \right] \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right) + \left(\varphi (T) - \theta \right) \right]^2 \\= \boldsymbol{\mathrm{E}} \left[ \left( \phi (\boldsymbol{X}) - \varphi (T) \right)^2 \right] + \boldsymbol{\mathrm{E}} \left[ \left(\varphi (T) - \theta \right)^2 \right] \\ \geq \boldsymbol{\mathrm{V}} \left[ \varphi (T) \right]

十分統計量に依存しない不偏推定量を用いた場合、下から2段目の式の第一項の分、推定量が悪くなり得ることを示しています。

まとめるとブラックウェルーラオ（Blackwell-Rao）の定理は以下のようになります[1]。

標本 \boldsymbol{X} = (X_1, \ldots, X_n)'の各要素が互いに独立同分布に従うとき（independent and identically distributed, i.i.d.と書くことも多い）、 T(\boldsymbol{X})を十分統計量、統計量 \phi (\boldsymbol{X})は母数 \thetaの不偏推定量とする。このとき

  1. \varphi (t) = \boldsymbol{\mathrm{E}} \left[ \phi (\boldsymbol{X}) | T=t \right] は不変推定量
  2. \boldsymbol{\mathrm{V}} \left[ \varphi (T) | \theta \right] \leq \boldsymbol{\mathrm{V}} \left[\phi (\boldsymbol{X}) \right] となる

参考書籍

  1. 応用を目指す数理統計学（国友 直人）