読者です 読者をやめる 読者になる 読者になる

カイ二乗検定

ピアソンのカイ二乗検定(Pearson's chi-square test)
カイ二乗検定のうち最も広く用いられている。
帰無仮説「観察された事象の相対的分布がある頻度分布に従う」という仮説を検定するもの。
観測値 O、帰無仮説から導き出される期待値 Eを用いて次式であらわされる。

 \chi^2 = \sum \frac{(O-E)^2}{E}

分散共分散行列

分散
まずは分散とは何かということから。
分散とは「各データが平均値からどれだけ離れているか」という、データのばらつき具合を表す。
具体的に、分散は「(各データの平均値からの距離)の二乗の平均」。
分散の平方根を取ったものが標準偏差と呼ばれるものである。
標準偏差 \sigmaで表すと、

 \sigma^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \overline{x})^2

 = \frac{1}{n} (X - \overline{X})^T (X - \overline{X})

これはある一つの事柄に対するデータがどれほど分散しているかを見るものである。

一方で、異なる事柄において相関があるのかを知りたいときに登場するのが共分散です。

共分散
分散共分散行列は次の式で表される。

 S_{XY} = \frac{1}{n} (X-\overline{X})^T (Y-\overline{Y})

この式が示すところを説明する。 たとえば、 Xが平均以上のとき、 Yも平均以上となる場合は、 X, Yの共分散値は正となり、
逆に Xが平均以上のとき、 Yは平均以下となる場合は、 X, Yの共分散値は負となる。

またそれぞれのデータに対して、「データから平均値を引き、標準偏差で割る」(基準化)と、分散共分散行列の対角成分が全て1になる。
つまり、元のデータに固有の平均値や標準偏差の大きさに影響されなくなる。
このようにして得られた行列は「相関行列」と呼ばれる。

分散共分散行列は、各成分が独立である(相関がない)場合は対角行列になる。 また分散共分散行列は半正定値である。 証明は以下の通り。

 S_{ij} = \frac{1}{N} \sum_{k=1}^{N} (x_i^{(k)} - m_i)(x_j^{(k)} - m_j) = \frac{1}{N} \sum_{k=1}^{N} z_i^{(k)} z_j^{(k)}
 S = \frac{1}{N} Z Z^T

任意のベクトル uに対して、 y=Z^{T} uとおくと

 u^T S u = \frac{1}{N} u^T Z Z^T u = \frac{1}{N} y^T y \geqq 0

ゆえに半正定値。

主成分分析(PCA)

今回は主成分分析(Principal Component Analysis)について勉強します。
これは統計データから互いに無関係の成分を取り出して、観測地をそれらの成分の線形結合で説明することを示す。

主成分分析の定式化
 n次元空間における主成分の軸を w_1, w_2, \cdots, w_nとする。
 xの各軸への写像 yは次式で与えられる。

 y_{1} = w_{11} x_{1} + w_{12} x_{2} + \cdots + w_{1n} x_{n} = w_1^{T} \vec{x}

 y_{2} = w_{21} x_{1} + w_{22} x_{2} + \cdots + w_{2n} x_{n} = w_2^{T} \vec{x}

 \vdots
 y_{n} = w_{n1} x_{1} + w_{n2} x_{2} + \cdots + w_{nn} x_{n} = w_n^{T} \vec{x}

ただし、

 w_{i}^{T} w_j = \delta_{ij}

ここでPCAとは、分散が最大となる方向ベクトル w_1, w_2, \cdots, w_nを求める問題である。

 y=W^T x

最小二乗法による主成分分析
平均を引いたサンプル点 x_j(j=1, 2, \cdots, N)とする。
 x_jの部分空間 w_1, w_2, \cdots, w_m (m \leq n)への写像 \hat{x}_j

 \hat{x}_j = y_1 w_1 + y_2 w_2 + \cdots + y_m w_m = \sum_{i=1}^{m} w_{i}^{T} x_j w_i

で与えられる。ここで、全てのサンプル点に対して、サンプル点とその写像との距離関数を次式で定義する。

 E(w_i) = \sum_{j=1}^{N} || x_j - \hat{x}_j ||
 = \sum_{j=1}^{N} || x_j - \sum_{i=1}^{m} w_i^T x_j w_i ||
 = \sum_{j=1}^{N}||x_j||^2 - 2 \sum_{j=1}^{N} \sum_{i=1}^{m} w_i^T x_j x_j^T w_i + \sum_{j=1}^{N} \sum_{i=1}^{m} w_i^T x_j x_j^T w_i
 = \sum_{j=1}^{N}||x_j||^2 - \sum_{j=1}^{N} \sum_{i=1}^{m} w_i^T x_j x_j^T w_i

ここで

 S = \sum_{j=1}^{N}x_j x_j^T

は共分散行列である。 ここで w_iの関数でない部分を取り除くと、 E(w_i)の最小化と

 \sum_{i=1}^{m} w_i^T S w_i

の最大化は同等である。したがってPCAは以下の拘束条件付き最適化問題となる。

 max \sum_{i=1}^{m} w_i^T S w_i

subject

 ||w_i||^2 = 1 (i=1, 2, \cdots, m)

これをLagrangeの未定乗数法を用いて解く。

 L(w_1, w_2, \cdots, w_m) = \sum_{i=1}^{m} w_i^T S w_i - \sum_{i=1}^{m} \lambda_i (w_i^T w_i - 1)

 w_iについて偏微分を取ると

 \frac{1}{2} \frac{\partial L}{\partial w_1} = S w_i - \lambda_i w_i = 0
 S w_i = \lambda_i w_i

以上から w_1, w_2, \cdots, w_mは共分散行列 S固有ベクトルとなる。 ここで yの共分散行列 S_y

 S_y = \sum_{j}^{N} y y^T = W^T S_x W = \begin{bmatrix}\lambda_1 & & & \\ & \lambda_2 & & \\ & & \ddots & \\ & & & \lambda_m \end{bmatrix}

したがって、共分散行列 S_x固有値問題を解くことに匹敵する。

ちなみに主成分分析はデータが正規分布に従っていることを仮定している。 そうでない場合を扱う場合、独立成分分析が主流となっている。

CommandLineParser[OpenCV]

以前にコマンドラインについての記事を書いたことがありました。
今回はOpenCVの関数を用いて、コマンドラインを読み込ませるプログラムを書いてみようと思います。

CommandLineParser
CommandLineParser、、、。
なんだチミはってか。そうです。こいつがコマンドラインを読み込むプルグラムを簡単に書く関数です。

CommandLineParser(int argc, char* argv[], const String &keys)

関数本体はこんな感じ。 argcとargvはmain分の引数と同じものをこの関数に渡し、この引数についての情報をkeysというやつに渡すということです。

Keysの書き方
Keysは3つの要素で構成されています。

  1. オプションの名前
  2. デフォルト値
  3. ヘルプメッセージ

OpenCVの公式documentの例を下に書きました。

const String keys = 
    "{help h usage ? |      | print this message   }"
    "{@image1        |      | image1 for compare   }"
    "{@image2        |<none>| image2 for compare   }"
    "{@repeat        |1     | number               }"
    "{path           |.     | path to file         }"
    "{fps            |-1.0  | fps for output video }"
    "{N count        |100   | count of objects     }"
    "{ts timestamp   |      | use time stamp       }"
    ;

このように「|」で区切って一つのオプションを書きます。
また一つのオプションに複数名前を付けたい場合はスペースで区切ります(例:help h usage ?)。
引数を位置指定するには前に@を付けます。
またデフォルト値のないオプションについてはhas()で確認できます。

CommandLineParser parser(argc, argv, keys);
if (parser.has("help"))
{
   parser.printMessage();
   return 0;
}

またStringの引数は次のように取得します。

std::string image1 = parser.get<std::string>("@image1");

また追加事項があれば随時更新します。 それでは。

上三角行列と下三角行列の性質

今回はこちらの記事こちらの文書を参考にします。

定義
 n \times nの正方行列 A=(a_{ij})について
 i > jならば a_{ij} = 0」を満たす行列を上三角行列(upper triangular matrix)
 j > iならば a_{ij} = 0」を満たす行列を下三角行列(lower triangular matrix)
という。

行列式の定義
三角行列の行列式を考える前に、一般の行列における行列式をおさらいしましょう。

 det(A) = \sum_{\sigma \in S_{n}} sgn(\sigma) \prod_{i=1}^{n} a_{n \sigma (i)}

 \sum_{\sigma \in S_n} sgn(\sigma) a_{1 \sigma (1)} a_{2 \sigma (2)} \cdots a_{n \sigma (n)}

  •  \sigmaは1から nの置換を表す。
  •  \sum_{\sigma \in S_n} n次の全ての置換に関して和を取ることを表す。
  •  sgn( \sigma )は置換の符号を表している。奇置換なら -1、偶置換なら 1をとる。

三角行列の行列式
上記した行列式の定義を利用すれば、逆行列を求めるのは容易です。
恒等置換以外では必ず 0が入るので、
 det(U) = \prod_{i=1}^{n} u_{ii}

三角行列の逆行列
ここでは上三角行列の逆行列が上三角行列であることを証明します。
まずは次の補題から証明します。
[補題1]

 A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\  & a_{22} & \cdots & a_{2n} \\ & & \ddots & \vdots \\ 0 & & & a_{nn} \end{bmatrix}

 n次の上三角行列、 a_{ii} \neq 0とすると、 AB=I_nを満たす上三角行列 B b_{ii} \neq 0が存在する。
[証明]
数学的帰納法により証明する。
 n=1のときは自明。
 n-1のときにこの補題が正しいとき、 n-1次の上三角行列 A' A' B' = I_nを満たす n-1次の上三角行列 B'が存在する。
ここで、

 B= \begin{bmatrix} b_{11} & b_{12} & \cdots & b_{1n-1} & x_{1} \\ & b_{22} & \cdots & b_{2n-1} & x_{2} \\ & & \ddots & \vdots & \vdots \\ & & & b_{n-1n-1} & x_{n-1} \\ 0 & & & 0 & x_{n} \end{bmatrix}


とおいて、 AB=I_nとなるような x_1, x_2, \cdots, x_nを求める。

 \left\{ \begin{array}{r} a_{11} x_{1} + a_{12} x_{2} + \cdots + a_{1n} x_{n} = 0 \\ a_{22} x_{2} + \cdots + a_{2n} x_{n} = 0 \\ \vdots \\ a_{nn} x_{n} = 1 \end{array} \right.

 a_{nn} \neq 0なので、

 x_n = a_{nn}^{-1}

これより

 x_{n-1} = - a_{n-1 n} a_{n n}^{-1} a_{n-1 n-1}^{-1}

となる。
以上の操作を繰り返すことにより、 x_1, x_2, \cdots, x_nを求めることができる。

[補題2]

 A= \begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\  & a_{22} & \cdots & a_{2n} \\ & & \ddots & \vdots \\ 0 & & & a_{nn} \end{bmatrix}

に対して、 n次行列 B=(b_{ij}) AB=I_nを満たすものが存在すれば、 a_{ii} \neq 0である。
[証明]
数学的帰納法で証明する。 n=1のときは自明。
 n-1のときにこの補題が成立するとき、

 AB_{n j} = \sum_{k = 1}^{n} a_{nk} b{kj} = a_{nn} b_{nj} = \delta_{n j}

よって、

 a_{nn} \neq 0

 i = 1, \cdots, n-1に対して

 a_{nj} = 0

が成立。
これを帰納法の仮定より、

 a_{ii} \neq 0

ちなみに余因子を使っても証明できます。
余因子ってなんだっけって人は下を参考にしてください。

余因子による逆行列
 A逆行列 ij成分は

 \Delta_{ij} / det(A)

 \Delta_{ij} A j行目と i列目を取り除いた行列の行列式 (-1)^{i + j}倍したもの。

特異値分解

今回はこちらのスライドこちらのスライドこちらのスライドを参考(ほぼ和訳)して特異値分解について理解したいと思います。

Introduction
$n$ 次元ベクトルの解を求めたいときに、誤差を含んだデータ $m(m>n)$ 個から解を得たいときのことを考えます。 画像認識(カメラキャリブレーションなど)では、線形方程式に対して解を計算する際にいくつかの方法がありますが、誤差を含むデータから解を得るにはロバスト推定をしたり、今回説明するように二乗誤差を最小化するようなことを考えます。

Motivation
線形方程式 $Ax=b$ に構成されるシステムを考えます。
ここで、 $b$ が $b+ \delta b$ といった具合に誤差が乗っているとき、解 $x$ には $A^{-1} \delta b$ の誤差が乗っかります。
そこで、 $|| \delta b ||$ によってどのような誤差が発生するかを以下の観点考えたい。 ( $ b $ はベクトルなので、 $\delta b$ ももちろんベクトルです。)
1. $b$ がどの方向に行けば誤差が大きくなるのか
2. 条件数 $cond(A)(= ||A|| \cdot ||A^{-1}||)$ が大きいとき、どのようにして正確な解を得ればよいか
上記のような問題について考えるときは特異値分解(Singular Value Decomposition)が有効である。

直行行列(Orthogonal Matrix)
 m \times n行列 U U^{T} U=Iを満たすとき、 $U$ は直行行列である。 つまり、 $U$ の列ベクトルは互いに直行している。

行列で見る特異値分解
$m \times n$ 行列 $A(m>n)$ のとき

  1.  m \times n 直行行列 $U$
     Uの列成分は AA^{T}固有ベクトルである。
     AA^{T} = USV^{T}VSU^{T} = US^{2}U^{T}
  2.  n \times n 対角行列  S
       S=diag (\sigma_1, \sigma_2, \cdots, \sigma_{n})とする。    \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_nのとき、 \sigma_1, \sigma_2, \cdots, \sigma_nを特異値と呼び、行列 A^{T}A固有値である。
  3. $n \times n$ 直行行列 $V$
     Vの列成分は A^{T}A固有ベクトルである。
     A^{T}A = VSU^{T}USDV^{T} = VS^{2}V^{T}

によって、
 A=USV^{T}
と分解することを $A$ を特異値分解するという。

 A= \sum_{j=0}^{n} \sigma_{j} u_{j} v_{j}^{T}

であるから、

 A v_{j} = \sigma_{j} u_{j}

となる。
この式の意味するところは、「 Vで張られる空間の基底 v_j Uで張られる空間の基底 u_jに射影される」ということである。
また、 \Lambda = S^{2}とすると、

 A A^{T}=U \Lambda U^{T}

であり、

 (A A^{T}) U = U \Lambda

 (A A^{T}) u_{j} = \lambda_{j} u_{j}

つまり、 Uの列成分は AA^{T}固有ベクトルであることが確認できる。
同様に

 A^{T} A=V \Lambda V^{T}

であり、

 (A^{T} A) V = V \Lambda

 (A^{T} A) v_{j} = \lambda_{j} v_{j}

つまり、  Vの列成分は A^{T}A固有ベクトルであることが確認できる。

固有値と特異値
一般に固有値と特異値は同じ値を取らない。
一般に、行列 Aに対して、 A^{H} Aはエルミート行列であり、非負の固有値を持つ。(ただし A^{H} Aのエルミート共役とする。)
特異値は A^{H} A固有値平方根である。

特異値分解による逆行列の計算
 n \times n行列 A A=USV^{T}特異値分解されるとき、正則である A逆行列
 A^{-1} = V (diag (\sigma_1^{-1}, \sigma_2^{-1}, \cdots, \sigma_n^{-1})) U^{T}
となる。
 \sigma_iがゼロに近いとき、丸め誤差が生じる。
ゼロに近い \sigma_iが多いほど、 Aは特異である(?)。  Aが特異であるまたは不良条件(ill-conditioned)であるとき、 A逆行列は、

 A^{-1} = V D^{-1}_{0} U^{T}

 D_{0}^{-1} = \left\{ \begin{array}{ll}
 1 / \sigma_i & (\sigma_i > t) \\
 0 & (otherwise) \end{array} \right.

ただし、 tは小さい閾値

The condition of a matrix
線形方程式 Ax=bにおいて
 bが微小変化したときに、解 xが比較的大きく変動するおき、不良条件(ill-conditioned)という。
最大特異値と最小特異値の比 \sigma_1 / \sigma_nが大きいほど、特異である。

ゼロ空間(null space)
核空間(kernel space)ともいう。
 n \times n行列 Aに対して、 Ax=0を満たすベクトル xの集合をゼロ空間(null space)という。

行列の値域(range)
行列 Aの値域(range)とは Ax=bの解 xが存在するベクトル bの集合のことである。
値域は線形ベクトル空間であり、次元は Aのランクと一致する。
つまり、ゼロでない特異値の数とランクが一致する。 もし Aが特異ならば、 Aのランクは nより小さい。
 n=Rank(A) + Nullity(A)

特異値分解と値域、ゼロ空間
特異値分解は値域とゼロ空間における直行基底を構成する。
 Aの非負の特異値に対応する、 Uの列ベクトルは、 Aの値域の正規直行基底である。
 Aのゼロの特異値に対応する、 Vの列ベクトルは、  Aのゼロ空間の直行基底を形成している。

線形方程式を特異値分解で解く
方程式 Ax=b(b \neq 0)であり、 b Aの値域に含まれるとき、 この連立方程式は解を持つ。
しかし、解が無数に存在することがある。(なぜならば、もし xが解で、 yがゼロ空間上のベクトルであるとき、 x+yも同様に解であるから。)
しかし、なんらかの値が必要になった場合、 ||x||^2が最も小さい解を解とする。
このとき解は
 x = V (diag (\sigma_1^{-1}, \sigma_2^{-1}, \cdots, \sigma_n^{-1})) (U^{T} b) = V S^{-1} U^{T} b
となる。ただし、 \sigma_{j}= 0のとき、 \sigma_j^{-1} = 0とする。

一方、 bが値域に含まれないとき、 Ax=bを満たすベクトル xは存在しなく、解は先ほどと同様に
 x = V (diag (\sigma_1^{-1}, \sigma_2^{-1}, \cdots, \sigma_n^{-1})) (U^{T} b) = V S^{-1} U^{T} b
で与えられる。
この解は ||Ax-b||が最小となるような解となっている。

なぜこれらのように bが値域に含まれても、含まれなくて同じように解が出せるかを証明する。
 bは値域に含まれるとき(つまり式の本数が足りないとき)
 Ax=bを拘束条件として、 ||x||を最小化させる。
ラグランジュの未定乗数法を用いて解く。

 L(x) = 1/2 ||x||^2 - (\lambda, Ax-b)

 \lambda = (\lambda_1, \lambda_2, \cdots, \lambda_m)

よって、

 x - A^{T} \lambda = 0

これによって、 Ax=b xに代入して、

 A A^{T} \lambda = b

これを解いて \lambdaを求めると、(つまり、 A A^{T}逆行列を持つことを前提としている)
次式を得る。

 x = A^T (A A^T)^{-1} b

次にこれが特異値分解による解と等しいことを示す。

 A A^{T} = U \Lambda U^{T}

これより

 A^{T} (A A^{T})^{-1} = V S U^{T} U \Lambda^{-1} U^{T}

 A^{T} (A A^{T})^{-1} = V S^{-1} U^{T}

となる。
次に、 bが値域に含まれないとき(つまり式の本数が多い時)

 1/2 ||Ax-b||^{2}

を最小化する。
これより、

 A^{T} A x - A^{T} b = 0

よって( A^{T} A逆行列を持つことを前提としている)

 x = (A^{T} A)^{-1} A^{T} b

となる。これが、特異値分解による解と等しいことを示す。

 (A^{T} A)^{-1} A^{T} = V \Lambda^{-1} V^{T} V S U^{T}
 A^{T} (A A^{T})^{-1} = V S^{-1} U^{T}

となる。

同次線形方程式を特異値分解で解く
同次線形方程式 Ax=0を考える。 これはゼロ空間上の任意のベクトルに等しい。
つまり、特異値がゼロに対応する Vの列ベクトルが作る空間が解となる。
(別の解釈?)←これが一番しっくりくる
最小ノルムとなる解 x=0は自明解。
そこで、 x=1という拘束条件を加える。
この拘束条件を満たす最適解を探索する問題に帰着。

  min_{||x||=1} ||Ax||^2

これをラグランジュの未定乗数法を使って解く。

 L(x) = x^{T} A^{T} A x  - \lambda (x^{T} x -1)

よって、

 A^{T} A x - \lambda x = 0

したがって、 \lambda A^{T} A固有値であり、そのときの固有ベクトル

 x= e_{\lambda}

とする。
ここで、 \lambda=0のとき L(e_{\lambda})を最小化できるので、 x=e_0が解となる。

(また別の解釈?)

 rank(A)=n-k (m \geq n-k, \sigma_{n-k+1}=\cdots=\sigma_n=0)

のとき、解は

 x=a_1 v_{n-k+1} + a_2 v_{n-k-1} + \cdots + a_k v_n

with

 a^{2}_{1} + a^{2}_{2} + \cdots + a^{2}_{k} = 1

Atomいろいろ

かなり久しぶりの投稿になります。お休みしていてすみませんでした。新年から(とはいえ今日は大晦日ですが)また定期的に投稿していこうと思います。

今回はAtomというテキストエディタについてです。 "AtomGitHubの創業者Chris Wanstrath氏が「Web技術を用いて、Emacsのように自由にカスタマイズできる新世代のエディターを開発する」という思いから始まったオープンソースのエディターです。"(こちらを参考) 素人の僕でも使えてかつなんかカッコイイという素晴らしいエディタです。GitHubとの連携もよくて、Markdownで書いたものをhtmlで保存できたりとか(これくらい普通か?)いい感じです!

インストール
何がともあれインストールから。

  1. 公式サイトからexeをインストール
  2. exeをインストールするとセットアップが始まります
  3. 言われるがままにボタンをポチポチ
  4. インストールが完了すると、atomが開きます

なんて簡単なんでしょう。簡単万歳。(こんなキャラでやってたっけ?)

Atomの利点
今回はとりあえず利点をいくつか挙げて詳細は後日追加していくことにします。
利点

  • リアルタイムでプレビューが見れる
  • markdownからHTMLへの変換可能
  • packageを利用することで、tex形式で数式を書くことができる
  • とにかく簡単

HTMLでMathjaxの設定
本記事での最終目標はここにしたい。覚え書です。無視してください。

<script type="text/javascript"
src="http://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML">
</script>