1.はじめに
もくじ
本稿は、「 カーネル関数の深層探訪：現代技術とその応用 」を参照しています。
1.1 カーネルトリック
- カーネルトリックの歴史は、機械学習の発展と密接に関連しています。
- この概念は、1960年代に統計学の分野で初めて導入されましたが、1990年代にサポートベクターマシン（SVM）とともに大きな注目を集めるようになりました。
- カーネルトリックの基本的なアイデアは、データをより高次元の特徴空間にマッピングすることで、線形分離が容易になるというものです。
- これは、非線形データセットに対して線形アルゴリズムを適用する際に特に有効です。
- この技術の発展には、VapnikとChervonenkisの研究が大きな役割を果たしました。
- 彼らは、SVMとともにカーネルメソッドの理論的基礎を築きました。
- カーネルトリックは、複雑な非線形問題を解決するための強力なツールとして、機械学習の分野で広く受け入れられています。
- 特に、画像認識、音声認識、自然言語処理など、多くの応用分野でその効果が実証されています。
- カーネルトリックの歴史を理解することは、現代の機械学習技術がどのように進化してきたかを理解する上で重要です。
- また、この技術が今後どのように発展し、新たな応用分野に適用される可能性があるかを考える上でも、その歴史的背景を知ることは有益です。
- 機械学習における線形と非線形の概念は、カーネルトリックを理解する上で不可欠です。
- 線形モデルは、入力特徴の線形組み合わせを使用して予測を行います。
- これは、データが直線や平面で分離可能な場合に適しています。
- しかし、多くの実世界の問題は非線形の性質を持ち、線形モデルでは適切に扱えません。
- 非線形モデルは、曲線や複雑な境界を用いてデータを分類します。
- これにより、より複雑なデータ関係を捉えることが可能になります。
- カーネルトリックは、この非線形性を扱うための鍵となります。
- カーネルトリックを使用することで、元の低次元空間の非線形関係を、高次元空間での線形関係に変換することができます。
- この変換により、線形アルゴリズムを用いて非線形問題を効率的に解決することが可能になります。
- このアプローチは、計算コストを抑えつつ、複雑なデータパターンを捉える能力を高めるため、機械学習において非常に価値が高いとされています。
- カーネル関数は、カーネルトリックの中核を成す要素であり、異なる種類のカーネルは異なる特性を持ちます。
- 代表的なカーネル関数には、線形カーネル、多項式カーネル、ガウシアン（RBF）カーネルなどがあります。
- 線形カーネルは最も単純で、元の特徴空間をそのまま使用します。
- これは計算コストが低いが、非線形問題には適していません。
- 多項式カーネルは、特徴の非線形組み合わせを可能にし、より複雑なデータパターンを捉えることができます。
- しかし、適切な多項式の次数を選ぶ必要があり、計算コストが高くなる可能性があります。
- ガウシアンカーネル（RBF）は、無限次元の特徴空間にデータをマッピングし、非常に複雑なデータ関係を捉えることができます。
- これは非線形問題に非常に効果的ですが、適切なパラメータの選択が重要です。
- カーネルの選択は、扱うデータの性質と特定の問題に依存します。適切なカーネル関数を選ぶことで、モデルの精度と効率を大幅に向上させることができます。
- したがって、異なるカーネル関数の理解と適用は、機械学習において重要なスキルとなります。
1.2 特殊な行列たち
- 以下のように、「左上～右下方向」に同じ値が並んだ行列を巡回行列（循環行列，Circulant matrix）という。
- 具体例として3次元ベクトル「$\vec{b}=(b_0,b_1,b_2)$」の場合だと次の通りです。
- 巡回行列の行列式は、因数分解が出来るので面白いらしい。
- 上記の巡回行列の行列式は、行列式の定義から以下のようになる。
- ここで「$\omega$」は、1の3乗根である（$\omega^3=1$）。
$\Large \displaystyle C= \begin{pmatrix} b_0 & b_1 \cdots & b_{N-2} & b_{N-1}\\ b_{N-1} & b_0 \cdots & b_{N-3} & b_{N-2}\\ \vdots & \vdots \ddots & \vdots & \vdots\\ b_2 & b_3 \cdots & b_0 & b_1\\ b_1 & b_2 \cdots & b_{N-1} & b_0 \end{pmatrix} $
$\Large \displaystyle C= \begin{pmatrix} b_0 & b_1 & b_2\\ b_2 & b_0 & b_1\\ b_1 & b_2 & b_0 \end{pmatrix} $
$\Large \displaystyle b_0^3 + b_1^3 + b_2^3 - 3b_0 b_1 b_2$
$\Large \displaystyle =(b_0 + b_1 + b_2)(b_0 + b_1\omega + b_2\omega^2)(b_0 + b_1\omega^2 + b_2\omega)$
- 置換行列とは、各行各列に一つだけ1の要素を持ち、それ以外全部ゼロの正方行列です。
- 具体例として4次元のとある置換行列は以下の通りです。
- （同じ4次元でもどこに1を当てはめるかで様々なパターンが存在します。）
- 置換行列は、掛けられたベクトルの要素を単純に入れ替えるだけの行列です。
- 試しに上記の具体例で示した置換行列を、4次元ベクトル「$\vec{x}=(x_0,x_1,x_2,x_3)$」にかけてみます。
- 単位行列「$I$」も掛けられたベクトルの各要素が元に戻ってくるとみなせるので、置換行列の特殊な場合です。
$\Large \displaystyle P= \begin{pmatrix} 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 \end{pmatrix} $
$\Large \displaystyle P\vec{x}= \begin{pmatrix} 0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 0 & 1\\ 0 & 0 & 1 & 0 \end{pmatrix} \begin{pmatrix} x_0\\ x_1\\ x_2\\ x_3 \end{pmatrix} = \begin{pmatrix} x_1\\ x_0\\ x_3\\ x_2 \end{pmatrix} $
- 巡回行列かつ置換行列のことを私が勝手に巡回置換行列と呼んでいます。
- つまり、要素のうちどれか一つだけ1、それ以外は全部ゼロのN次元の単位ベクトル「$\vec{e}=(0,1,0,\cdots,0,0)$」を巡回させながら構成した置換行列です。
- FIRフィルタによる畳み込み積分を行列で表現すると巡回行列になり、その巡回行列をさらに分解するとこの巡回置換行列で構成されていることがわかりました。
- 今回の話はここから始まります。
$\Large \displaystyle P= \begin{pmatrix} 0 & 1 & 0 & \cdots & 0 & 0\\ 1 & 0 & 0 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots\\ 0 & 0 & 0 & \cdots & 1 & 0\\ 0 & 0 & 0 & \cdots & 0 & 1\\ 1 & 0 & 0 & \cdots & 0 & 0 \end{pmatrix} $
- 巡回行列「$C$」の1行目を「$(c_0,c_1,\cdots,c_{n-1})$」とする。
- 巡回行列「$C$」の固有ベクトル「$\vec{x_k}$」は、以下のようになる。
- 巡回行列「$C$」の固有値「$\lambda_k$」は、以下のようになる。ただし、「$k=0,1,\cdots,n-1$」とする。
- 巡回行列「$C$」の固有値「$\lambda_k$」と固有ベクトル「$\vec{x}_k$」の関係は、以下のようになる。
- n=3の場合は以下のようになる。
- k=0:
- k=1:
- k=2:
- 従って、巡回行列「$C$」の行列式は、「$C$」の1行目を「$(c_0,c_1,\cdots,c_{n-1})$」とすると、以下のようになる。
- ここで、「$\zeta_{k}$」は、1のn乗根である（$\zeta_{k}^n=1$）。
$\Large \displaystyle \vec{x_k}=(1,\lambda_k,\lambda_k^2,\cdots,\lambda_k^{n-1})^T$
$\Large \displaystyle \lambda_k=\sum_{l=0}^{n-1}\zeta_{k}^{l}c_{l}\ \ \ (定理2)$
$\Large \displaystyle C \vec{x}_k=\lambda_k \vec{x}_k$
$\Large \displaystyle \begin{pmatrix} c_0 & c_1 & c_2\\ c_2 & c_0 & c_1\\ c_1 & c_2 & c_0 \end{pmatrix} \begin{pmatrix} 1\\ 1\\ 1 \end{pmatrix} = (c_0+c_1+c_2) \begin{pmatrix} 1\\ 1\\ 1 \end{pmatrix} $
$\Large \displaystyle \lambda_0=c_0+c_1+c_2$
$\Large \displaystyle \begin{pmatrix} c_0 & c_1 & c_2\\ c_2 & c_0 & c_1\\ c_1 & c_2 & c_0 \end{pmatrix} \begin{pmatrix} 1\\ \omega\\ \omega^2 \end{pmatrix} = (c_0+\omega c_1+\omega^2 c_2) \begin{pmatrix} 1\\ \omega\\ \omega^2 \end{pmatrix} $
$\Large \displaystyle \lambda_1=c_0+w c_1+w^2 c_2$
$\Large \displaystyle \begin{pmatrix} c_0 & c_1 & c_2\\ c_2 & c_0 & c_1\\ c_1 & c_2 & c_0 \end{pmatrix} \begin{pmatrix} 1\\ \omega^2\\ \omega \end{pmatrix} = (c_0+\omega^2 c_1+\omega c_2) \begin{pmatrix} 1\\ \omega^2\\ \omega \end{pmatrix} $
$\Large \displaystyle \lambda_2=c_0+w^2 c_1+w c_2$
$\Large \displaystyle det C=\prod_{k=0}^{n-1}\lambda_k$
$\Large \displaystyle =\prod_{k=0}^{n-1}(\sum_{l=0}^{n-1}\zeta_{k}^{l}c_{l})$
$\Large \displaystyle \zeta_{k}=exp(\frac{2\pi kj}{n})$
$\Large \displaystyle \zeta_{k}^n=exp(\frac{2\pi kj+2\pi kj+,\cdots,+2\pi kj}{n})$
$\Large \displaystyle =exp(2\pi kj)=1$
- 「$N\times N$」の巡回置換行列「$P$」を単位行列「$I$」に「$N$」回かけると、また元の単位行列「$I$」に戻ってくる事がわかります。
- 4次元の巡回置換行列で試してみます。
1.3 合成積
- 二つの関数「$f(x)$」と「$g(x)$」から以下のような操作をして新しい関数「$h(x)$」を作る。
- この「$h(x)$」を「$f(x)$」と「$g(x)$」の合成積（畳み込み）と呼び、以下のように表記できる。
$\Large \displaystyle h(x)=\int_{-\infty}^{\infty}f(t)g(x-t)dt$
$\Large \displaystyle h(x)=f(x)*g(x),\ \ h=f*g$
- 二つの数列「$\{a_n\}$」と「$\{b_n\}$」から以下のような操作をして新しい数列「$\{c_n\}$」を作る。
- この「$\{c_n\}$」を「$\{a_n\}$」と「$\{b_n\}$」の合成積（畳み込み）と呼び、以下のように表記できる。
- 例えば「$n=6$」の時、合成積「$c_6$」は、以下のようになる。
$\Large \displaystyle c_n=\sum_{t=0}^{n}a_t b_{n-t}$
$\Large \displaystyle c_n=a_n*b_n,\ \ c=a*b$
$\Large \displaystyle c_6=a_0 b_6 + a_1 b_5 + a_2 b_4 + a_3 b_3 + a_4 b_2 + a_5 b_1 + a_6 b_0$
- 回転因子「$w$」を以下のように定義する。
- 「$w^a$」は、以下の例のように単位円（半径1の円）に現れる。
$\Large \displaystyle w=e^{-\frac{2\pi}{N}i}$
- 回転因子「$w$」を使うと、離散フーリエ変換は以下のようになる。
- 回転因子「$w$」を使って、離散フーリエ変換「$F(k)$」を行列で表す。
- なので、行列で表すと以下のようになる。
- DFT行列「$W_N$」を使用する、以下のように表す事ができる。
- DFT行列「$W_N$」を使用する、以下のようにフーリエ変換を線形変換で表す事ができる。
- 例えば、「$n=3$」の場合、DFT行列「$W_N$」は、次のようになる。
$\Large \displaystyle F(k)=\sum_{n=0}^{N-1}f(n)e^{-i\frac{2\pi k}{N}n}$
$\Large \displaystyle =\sum_{n=0}^{N-1}f(n)e^{(-\frac{2\pi}{N})kn}$
$\Large \displaystyle =\sum_{n=0}^{N-1}f(n)w^{kn}$
$\Large \displaystyle F(0)=w^{0-0}f(0)+w^{0-1}f(1)+,\cdots,+w^{0-(N-1)}f(N-1)$
$\Large \displaystyle F(1)=w^{1-0}f(0)+w^{1-1}f(1)+,\cdots,+w^{1-(N-1)}f(N-1)$
$\Large \displaystyle F(2)=w^{2-0}f(0)+w^{2-1}f(1)+,\cdots,+w^{2-(N-1)}f(N-1)$
$\Large \displaystyle F(N-1)=w^{(N-1)-0}f(0)+w^{(N-1)-1}f(1)+,\cdots,+w^{(N-1)-(N-1)}f(N-1)$
$\Large \displaystyle \begin{pmatrix} F(0) \\ F(1) \\ F(2) \\ \vdots \\ F(N-1) \end{pmatrix} = \begin{pmatrix} w^{0-0} & w^{0-1} & w^{0-2} & \cdots & w^{0-(N-1)} \\ w^{1-0} & w^{1-1} & w^{1-2} & \cdots & w^{1-(N-1)} \\ w^{2-0} & w^{2-1} & w^{2-2} & \cdots & w^{2-(N-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ w^{(N-1)-0} & w^{(N-1)-1} & w^{(N-1)-2} & \cdots & w^{(N-1)-(N-1)} \end{pmatrix} \begin{pmatrix} f(0) \\ f(1) \\ f(2) \\ \vdots \\ f(N-1) \end{pmatrix} $
$\Large \displaystyle \begin{pmatrix} F(0) \\ F(1) \\ F(2) \\ \vdots \\ F(N-1) \end{pmatrix} = W_N \begin{pmatrix} f(0) \\ f(1) \\ f(2) \\ \vdots \\ f(N-1) \end{pmatrix} $
$\Large \displaystyle \vec{y}=W_N \vec{x}$
$\Large \displaystyle W_N= \begin{pmatrix} 1 & 1 & 1 \\ 1 & w & w^2 \\ 1 & w^2 & w^4 \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ 1 & w & w^2 \\ 1 & w^2 & w \end{pmatrix} $
- 「定理2」では、巡回行列の固有ベクトルを求めたが、これらを並べた行列は、離散フーリエ変換を表す行列 「$W_N$」と一致する。
- 例えば、「$n=3$」の場合、3次元ベクトル 「$\vec{x}$」の離散フーリエ変換は、以下のように表せる。
- 巡回行列「$C$」の固有値「$\lambda_k$」と固有ベクトル「$\vec{x}_k$」の関係は、以下のようになる。
- 例えば、「$n=3$」の場合、に表せる。
$\Large \displaystyle W_N\vec{x}= \begin{pmatrix} 1 & 1 & 1\\ 1 & \omega & \omega^2\\ 1 & \omega^2 & \omega^4 \end{pmatrix} \vec{x} = \begin{pmatrix} 1 & 1 & 1\\ 1 & \omega & \omega^2\\ 1 & \omega^2 & \omega \end{pmatrix} \vec{x} $
$\Large \displaystyle C \vec{x}_k=\lambda_k \vec{x}_k$
$\Large \displaystyle \vec{x}_k=(1,\lambda_k,\lambda_k^2,\cdots,\lambda_k^{N-1})^T$
$\Large \displaystyle C \vec{x}_0=\lambda_0 \vec{x}_0$
$\Large \displaystyle C \vec{x}_1=\lambda_1 \vec{x}_1$
$\Large \displaystyle C \vec{x}_2=\lambda_2 \vec{x}_2$
以上
2.線形回帰モデル
もくじ
- 2.1 単回帰
- 2.2 重回帰とベクトル表現
- 2.3 線形回帰モデル
- 2.4 リッジ回帰
- 2.5 相関係数と回帰モデル
- 2.6 HSIC (Hilbert-Schmidt Independence Criterion)
2.1 単回帰
- ある入力「$x$」が与えられた時、出力「$y$」を予想する事を「回帰(regression)」と言う。
- 「ガウス過程」も「回帰(regression)」のための道具である。
- 最も単純な単回帰を考える。
- 単回帰では、以下のような一次式で、入力「$x$」から出力「$y$」の予測値「$\hat{y}$」を推定する。
- この時、パラメータ「$a$」と「$b$」をどうやって決めるかが単回帰の問題である。
- で、説明変数と被説明変数の関係性をモデル化する。、
- 予測値を「$\hat{y}=a+b x_n$」を書く事にすると、実際の「$y_n$」との誤差は、次のようになる。
- 上記の誤差を最小限にするために、以下の評価関数「$E$」を最小にするパラメータ「$a,b$」を求める。
- 上式が最小になるには、上式の「$a$」の偏微分と「$b$」の偏微分とが0になれば良い。
- だから、「$a$」の偏微分と「$b$」の偏微分は、以下のようになる。
- 従って、以下の式を得る。
- (1)式から、以下の「$a$」を得る。
- 上式の「$a$」を(2)式に代入する。
- 上式から、以下の「$b$」を得る。
- 従って、以下の「$a$」を得る。
- 例として「$\{(3,2),(2,4),(-1,1)\}$」を計算して見る。
- 従って、パラメータ「$a,b$」は、以下のようになる。
$\Large \displaystyle \hat{y}=a+bx\ \ (a,b\in R)$
$\Large \displaystyle y_n - \hat{y}_n=y_n - (a+bx_n)$
$\Large \displaystyle E=\sum_{n=1}^{N}(y_n - \hat{y}_n)^2=\sum_{n=1}^{N}(y_n - (a+b x_n))^2$
$\Large \displaystyle E=\sum_{n=1}^{N}(y_n - a - b x_n)^2$
$\Large \displaystyle =\sum_{n=1}^{N}(y_n^2 + a^2 + b^2 x_n^2 + 2a y_n -2a b x_n - 2b x_n y_n)$
$\Large \displaystyle \frac{\partial E}{\partial a}=\sum_{n=1}^{N}(2a - 2y_n + 2b x_n)=0$
$\Large \displaystyle \frac{\partial E}{\partial b}=\sum_{n=1}^{N}(2b x_n^2 + 2a x_n - 2x_n y_n)=0$
$\Large \displaystyle a N + b \sum_{n=1}^{N} x_n=\sum_{n=1}^{N} y_n\ \ (1)$
$\Large \displaystyle a \sum_{n=1}^{N}x_n + b \sum_{n=1}^{N} x_n^2=\sum_{n=1}^{N} x_n y_n\ \ (2)$
$\Large \displaystyle a= \frac{1}{N}(\sum_{n=1}^{N} y_n - b \sum_{n=1}^{N} x_n)=\bar{y}-b\bar{x}$
$\Large \displaystyle \sum_{n=1}^{N}x_n(\sum_{n=1}^{N} y_n - b \sum_{n=1}^{N} x_n) + N b \sum_{n=1}^{N} x_n^2=N\sum_{n=1}^{N} x_n y_n$
$\Large \displaystyle b (-(\sum_{n=1}^{N}x_n)^2 + N \sum_{n=1}^{N} x_n^2)=N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n$
$\Large \displaystyle b=\frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$
$\Large \displaystyle a=\frac{\sum_{n=1}^{N} y_n}{N} - \frac{\sum_{n=1}^{N} x_n}{N} \frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$
$\Large \displaystyle =\frac{\sum_{n=1}^{N} y_n}{N}\frac{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2} - \frac{\sum_{n=1}^{N} x_n}{N} \frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$
$\Large \displaystyle = \frac{\sum_{n=1}^{N} x_n^2\sum_{n=1}^{N} y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} x_n y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$
$\Large \displaystyle \sum_{n=1}^{N} x_n=3+2-1=4$
$\Large \displaystyle \sum_{n=1}^{N} y_n=2+4+1=7$
$\Large \displaystyle \sum_{n=1}^{N} x_n^2=9+4+1=14$
$\Large \displaystyle \sum_{n=1}^{N} x_n y_n=6+8-1=13$
$\Large \displaystyle a=\frac{14*7-4*13}{3*14-4*4}=\frac{98-52}{42-16}=\frac{46}{26}=1.77$
$\Large \displaystyle b=\frac{3*13-4*7}{3*14-4*4}=\frac{39-28}{42-16}=\frac{11}{26}=0.42$
- まず、2☓2行列Aの逆行列Bを求めます。
- 「$AB=I$」（単位行列）だから、
- 従って、以下の①～④の式を得る。、
- 式①☓d-式③☓bからeを得る。
- 式②☓d-式④☓bからfを得る。
- 式①☓c-式③☓aからgを得る。
- 式②☓c-式④☓aからhを得る。
- これを行列表示すると、以下のAの逆行列Bを得る。
$\Large \displaystyle A= \begin{pmatrix} a & b \\ c & d \end{pmatrix} $
$\Large \displaystyle B= \begin{pmatrix} e & f \\ g & h \end{pmatrix} $
$\Large \displaystyle \begin{pmatrix} a & b \\ c & d \end{pmatrix} \begin{pmatrix} e & f \\ g & h \end{pmatrix} = \begin{pmatrix} ae+bg & af+bh \\ ce+dg & cf+dh \end{pmatrix} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} $
$\Large \displaystyle ae+bg=1\ \ ①$
$\Large \displaystyle af+bh=0\ \ ②$
$\Large \displaystyle ce+dg=0\ \ ③$
$\Large \displaystyle cf+dh=1\ \ ④$
$\Large \displaystyle (ae+bg)d-(ce+dg)b=ade-bce=e(ad-bc)=d$
$\Large \displaystyle e=\frac{d}{ad-bc}$
$\Large \displaystyle (af+bh)d-(cf+dh)b=adf-bcf=f(ad-bc)=-b$
$\Large \displaystyle e=\frac{-b}{ad-bc}$
$\Large \displaystyle (ae+bg)c-(ce+dg)a=bcg-adg=-g(ad-bc)=c$
$\Large \displaystyle g=\frac{-c}{ad-bc}$
$\Large \displaystyle (af+bh)c-(cf+dh)a=bch-adh=-h(ad-bc)=-a$
$\Large \displaystyle h=\frac{a}{ad-bc}$
$\Large \displaystyle A^{-1}=B=\frac{1}{ad-bc} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} = \frac{1}{|A|} \begin{pmatrix} d & -b \\ -c & a \end{pmatrix} $
- ここでは、今後の見通しを考えて「(1)式と(2)式」を以下のような正規方程式で表現して見る。
- 「3.2☓2行列の逆行列」で求めた以下の式を利用して「$(a b)^T$」を求める。
- 「$(a b)^T$」を求める。
- となり、「2.単回帰」の偏微分で得られたパラメータ「a」と「b」の値と同じになる。
$\Large \displaystyle \begin{pmatrix} N & \sum_{n=1}^{N} x_n \\ \sum_{n=1}^{N} x_n & \sum_{n=1}^{N} x_n^2 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} = \begin{pmatrix} \sum_{n=1}^{N} y_n \\ \sum_{n=1}^{N} x_n y_n \end{pmatrix} $
$\Large \displaystyle \begin{pmatrix} A & B \\ C & D \end{pmatrix} ^{-1}= \frac{1}{AD-BC} \begin{pmatrix} D & -B \\ -C & A \end{pmatrix} $
$\Large \displaystyle \begin{pmatrix} a \\ b \end{pmatrix} = \frac{1}{N\sum_{n=1}^{N} x_n^2-(\sum_{n=1}^{N} x_n)^2} \begin{pmatrix} \sum_{n=1}^{N}x_n^2 & -\sum_{n=1}^{N} x_n\\ -\sum_{n=1}^{N}x_n & N \end{pmatrix} \begin{pmatrix} \sum_{n=1}^{N}y_n \\ \sum_{n=1}^{N}x_n y_n \end{pmatrix} $
$\Large \displaystyle = \frac{1}{N\sum_{n=1}^{N} x_n^2-(\sum_{n=1}^{N} x_n)^2} \begin{pmatrix} \sum_{n=1}^{N}x_n^2 \sum_{n=1}^{N}y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N}x_n y_n\\ N\sum_{n=1}^{N}x_n y_n - \sum_{n=1}^{N}x_n \sum_{n=1}^{N}y_n \end{pmatrix} $
2.2 重回帰とベクトル表現
- 重回帰では、入力「$x$」が一次元ではなく、次のようなD次元のベクトルになると考える。
- 重回帰では、入力「$x$」が一次元ではなく、次のようなD次元のベクトルになると考える。
- 単回帰の場合のように、二乗誤差を求めると以下のようになるが、このまま微分するのは見通しが悪い。
- そこで、以下のような係数ベクトル「$w$」を導入する。
- 更に、最初の係数が1になりような入力データベクトル「$w$」を導入する。
- 予測値「$\hat{y}$」は、以下のように「$w$」と「$x$」の内積で表現できる。
- 上式をN個のデータについて並べると以下のようになる。
- 上式をN個のデータについて並べると以下のようになる。
- 以下のように行列「$X$」と「$\hat{y}$」を定義すると、
- 「$\hat{y}$」は、以下のように書ける。この時、行列「$X$」を正則行列と呼ぶ。
$\Large \displaystyle X^T=(x_1,x_2,\cdots,x_D)$
$\Large \displaystyle \hat{y}=w_0+w_1 x_1+w_2 x_2\cdots+w_D x_D)$
$\Large \displaystyle (y_n-\hat{y}_n)^2=(y_n - (w_0+w_1 x_1+w_2 x_2\cdots+w_D x_D))$
$\Large \displaystyle w^T=(w_0,w_1,w_2,\cdots,w_D)$
$\Large \displaystyle x^T=(1,x_1,x_2,\cdots,x_D)$
$\Large \displaystyle \hat{y}=w_0+w_1 x_1+w_2 x_2\cdots+w_D x_D)$
$\Large \displaystyle =(w_0,w_1,w_2,\cdots,w_D) \begin{pmatrix} 1 \\ x_1 \\ x_2 \\ \vdots \\ x_D \end{pmatrix} =w^T x $
$\Large \displaystyle \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} w^T x_1 \\ w^T x_2 \\ \vdots \\ w^T x_N \end{pmatrix} = \begin{pmatrix} x_1^T \\ x_2^T \\ \vdots \\ x_N^T \end{pmatrix} w $
$\Large \displaystyle \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} = \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1D} \\ 1 & x_{21} & x_{22} & \cdots & x_{2D} \\ \vdots & & & & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{ND} \\ \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ w_2 \\ \vdots \\ w_D \end{pmatrix} $
$\Large \displaystyle X= \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1D} \\ 1 & x_{21} & x_{22} & \cdots & x_{2D} \\ \vdots & & & & \vdots \\ 1 & x_{N1} & x_{N2} & \cdots & x_{ND} \\ \end{pmatrix} $
$\Large \displaystyle \hat{y}= \begin{pmatrix} \hat{y}_1 \\ \hat{y}_2 \\ \vdots \\ \hat{y}_N \end{pmatrix} $
$\Large \displaystyle \hat{y}=Xw $
2.3 線形回帰モデル
2.4 リッジ回帰
2.5 相関係数と回帰モデル
- 回帰モデルと相関係数には、密接な関係がある。
- xの平均値「\bar{x}\」とyの平均値「\bar{y}\」は、以下のように定義される。
- すると単回帰モデルの直線「\y=a+bx$」の傾き「$b$」の最小二乗解は、以下のようになる。
- 一方、
- なので、傾き「$b$」は、
- とも書ける。
- ここで、「$\gamma_{xy}$」は、ピアソンの積率相関係数であり、 「$\sigma_{x}$」と「$\sigma_{y}$」は、それぞれ変数「$x$」「$y$」の標準偏差である。
- 「$\frac{\sigma_y}{\sigma_x}$」は、変数「$x$」「$y$」のスケールの比を表す定数であり、 変数「$x$」「$y$」の偏差が等しくなるように事前にデータを正規化しておけば、この値は1になる。
- 従って、傾き「$b$」は、実質的に積率相関係数「$\gamma_{xy}$」を表している。
$\Large \displaystyle \bar{x}=\frac{1}{N}\sum_{n=1}^{N}x_n$
$\Large \displaystyle \bar{y}=\frac{1}{N}\sum_{n=1}^{N}y_n$
$\Large \displaystyle b=\frac{N\sum_{n=1}^{N} x_n y_n - \sum_{n=1}^{N}x_n\sum_{n=1}^{N} y_n}{N \sum_{n=1}^{N} x_n^2 - (\sum_{n=1}^{N}x_n)^2}$
$\Large \displaystyle =\frac {\sum_{n=1}^{N} x_n y_n - N\sum_{n=1}^{N}\bar{x}\bar{y}} {\sum_{n=1}^{N} x_n^2 - N\bar{x}^2}$
$\Large \displaystyle \sum_{n=1}^{N}(x_n-\bar{x})(y_n-\bar{y}) = \sum_{n=1}^{N}(x_n y_n - x_n \bar{y}- \bar{x}y_n + \bar{x}\bar{y}) = \sum_{n=1}^{N}(x_n y_n) - N\bar{x}\bar{y} $
$\Large \displaystyle \sum_{n=1}^{N}(x_n-\bar{x})^2 = \sum_{n=1}^{N}(x_n^2 - 2\bar{x}x_n + \bar{x}^2) = \sum_{n=1}^{N}(x_n^2) - N\bar{x}^2 $
$\Large \displaystyle b= \frac{\sum_{n=1}^{N}(x_n-\bar{x})(y_n-\bar{y})} {\sum_{n=1}^{N}(x_n-\bar{x})^2} $
$\Large \displaystyle = \frac{\sum_{n=1}^{N}(x_n-\bar{x})(y_n-\bar{y})} {\sqrt{\sum_{n=1}^{N}(x_n-\bar{x})^2}\sqrt{\sum_{n=1}^{N}(y_n-\bar{y})^2}} \cdot \frac{\sqrt{\sum_{n=1}^{N}(y_n-\bar{y})^2}} {\sqrt{\sum_{n=1}^{N}(x_n-\bar{x})^2}} $
$\Large \displaystyle =\gamma_{xy}\cdot \frac{\sigma_y}{\sigma_x} $
2.6 HSIC (Hilbert-Schmidt Independence Criterion)
- HSIC (Hilbert-Schmidt Independence Criterion：Hilbert-Schmidt独立性基準)は、 2変数間の依存関係を評価する手法の一つであり、従来のピアソンの積率相関係数が変数間の線形な関係しか捉えることができない欠陥があり、 非線形な関係を取り扱う指標としてHSICを用いた検定が考案された。
- HSICは2変数「$X$」と「$Y$」の間の独立性を評価するノンパラメトリックな指標で、 より値が0に近いほど2変数は独立であると言うことができる。
- イメージとしては、複雑な非線形関係にある「$X$」と「$Y$」も、何かとっても都合の良い写像で飛ばした先では 単純な関係であり、独立性も簡単に評価できるんじゃないかという想いから来ているようだ。
- HSICでは、再生核ヒルベルト空間にデータを飛ばし、そこでcorss-covariance operatorのノルム(HS norm)を とった値を用い、permutation testで独立性を評価するらしい.
- HSICの定義とempirical estimatorを与えたことがGretton et al. (2005)の功績らしい。
以上
3.ガウス分布
もくじ
3.1 一次元正規分布
- 以下の確率密度関数「$p(x)$」で与えられる確率変数を正規分布と言う。
- 「$\mu$」は平均値であり、「$\sigma$」は標準偏差であり、「$\sigma^2$」は分散である。
- 上記のような正規分布を「$N(\mu,\sigma^2)$」と表記する。「$N(0,1)$」は「標準正規分布」と呼ばれる。
$\Large \displaystyle p(x)=\frac{1}{\sqrt{2\pi\sigma^2}}exp\{-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}\}$
- ①1次モーメント（平均）と2次モーメント（分散）を与える事により、それに従う確率法則を完全に確定できる。
- ②確率変数「$X$」が正規分布「$N(\mu,\sigma^2)$」に従う時、その線形変換「$Y=aX+b$」は「$N(a\mu+b,a^2\sigma^2)$」に従う。
- ③
- ④正規分布の確率密度関数をフーリエ変換すると、やはり正規分布の関数になる。
3.2 多次元正規分布
- ①期待値ベクトルと共分散行列を指定すれば分布が定まる。
- ②相関がないだけで独立だと断言できる。
- ③多次元正規分布を線形変換したら、また多次元正規分布になる。
3.3 共分散と相関係数
- 確率変数「$X$」に対する期待値「$\mu$」は、以下のように定義される。
- 確率変数「$X$」に対する分散「$V[X]$」は、以下のように定義される。
- 確率変数「$X$」に対する標準偏差「$\sigma$」は、以下のように定義される。
$\Large \displaystyle \mu=E[X]$
$\Large \displaystyle V[X]=E[(X-\mu)^2]$
$\Large \displaystyle \sigma=\sqrt{V[X]}$
- 確率変数「$X,Y$」に対する共分散「$Cov[X,Y]$」を導入する。
- 確率変数「$X,Y$」の期待値をそれぞれ「$\mu,\nu$」とする時、「$X,Y$」の共分散を以下のように定義する。
$\Large \displaystyle Cov[X,Y]=E[(X-\mu)(Y-\nu)]$
- 共分散の定義から、以下のような共分散の性質があります。
$\Large \displaystyle Cov[X,Y]=Cov[Y,X]$
$\Large \displaystyle Cov[X,X]=V[X]$
$\Large \displaystyle Cov[X+a,Y+b]=Cov[X,Y]$
$\Large \displaystyle Cov[aX,bY]=abCov[X,Y]$
- 上記の共分散の性質から、確率変数「$X,Y$」に定数「$a,b$」を掛けたら共分散が「$ab$」倍になると言う問題があった。
- そこで、「尺度」に惑わされない指標として以下の相関係数「$\rho_{XY}$」を導入する。
- つまり、以下のように確率変数「$X,Y$」を分散の平方根（標準偏差）で割って、正規化する。
- つまり、以下のように確率変数「$X,Y$」を分散の平方根（標準偏差）で割って、正規化したものの共分散が相関係数になる。
$\Large \displaystyle \rho_{XY}=\frac{Cov[X,Y]}{\sqrt{\mathstrut V[X]}\sqrt{\mathstrut V[Y]}}$
$\Large \displaystyle \hat{X}=\frac{X}{\sqrt{V[X]}}$
$\Large \displaystyle \hat{Y}=\frac{Y}{\sqrt{V[Y]}}$
$\Large \displaystyle Cov[\hat{X},\hat{Y}]=Cov[\frac{X}{\sqrt{V[X]}},\frac{Y}{\sqrt{V[Y]}}]=\frac{Cov[X,Y]}{\sqrt{\mathstrut V[X]}\sqrt{\mathstrut V[Y]}}=\rho_{XY}$
3.4 共分散行列
- 例えば、以下のように確率変数「$X_1,X_2,X_3$」の全ての組み合わせの共分散の一覧表を作り、行列化する。
- ここで、確率変数「$X_1,X_2,\cdots,X_n$」を縦に並べた縦ベクトル「$X$」を導入すると、その分散「$V[X]$」は、以下のように行列化された共分散一覧になる。
| $X_1$ | $X_2$ | $X_3$ | |
|---|---|---|---|
| $X_1$ | $Cov(X_1,X_1)$ | $Cov(X_1,X_2)$ | $Cov(X_1,X_3)$ |
| $X_2$ | $Cov(X_2,X_1)$ | $Cov(X_2,X_2)$ | $Cov(X_2,X_3)$ |
| $X_3$ | $Cov(X_3,X_1)$ | $Cov(X_3,X_2)$ | $Cov(X_3,X_3)$ |
$\Large \displaystyle V[X]= \begin{pmatrix} V[X_1] & Cov[X_1,X_2] & Cov[X_1,X_3] \\ Cov[X_2,X_1] & V[X_2] & Cov[X_2,X_3] \\ Cov[X_3,X_1] & Cov[X_3,X_2] & V[X_3] \\ \end{pmatrix} $
- 共分散行列を宣言しただけでは、あまりご利益はないので、確率変数「$X_1,X_1,\cdots,Xn$」を縦に並べた、 以下のような縦ベクトル「$X$」を導入する。
- 縦ベクトル「$X$」を使うと、共分散行列を以下のようにまとめて書ける。
- ただし、「$\mu$」は、以下のように定義される。
- つまり、縦ベクトル「$X$」の共分散行列（「$X$」の分散「$V[X]$」）は以下のようになる。
$\Large \displaystyle X= \begin{pmatrix} X_1 \\ X_2 \\ \vdots \\ X_n \end{pmatrix} $
$\Large \displaystyle V[X]=E[(X-\mu)(X-\mu)^T]$
$\Large \displaystyle \mu=E[X]$
$\Large \displaystyle V[X]= \begin{pmatrix} V[X_1] & Cov[X_1,X_2] & Cov[X_1,X_3] \\ Cov[X_2,X_1] & V[X_2] & Cov[X_2,X_3] \\ Cov[X_3,X_1] & Cov[X_3,X_2] & V[X_3] \\ \end{pmatrix} $
3.5 自己共分散の定義
$\Large \displaystyle \gamma_k =\frac{1}{T}\sum_{t=k+1}^{T}(y_t-\bar{y})(y_{t-k}-\bar{y})$
4.ガウス過程
もくじ
4.1 線形回帰モデルと次元の呪い
- 線形回帰モデルでは、基底関数による特徴ベクトル「$\phi(x)$」を導入する事で複雑な関数「$y$」が表現できる。
- まず、一次元の入力「$x$」について、以下のような特徴ベクトル「$\phi(x)$」と重みベクトル「$w$」を導入する。
- 例えば、以下のような一次元の入力「$x$」についての3次式の関数「$y$」は、 以下のように重みベクトル「$w$」と特徴ベクトル「$\phi(x)$」の内積で表す事ができる。
$\Large \displaystyle \phi(x)=(1,x,x^2,x^3)^T$
$\Large \displaystyle w=(w_0,w_1,w_2,w_3)^T$
$\Large \displaystyle y=w_0+w_1 x + w_2 x^2 + w_3 x^3$
$\Large \displaystyle =w^T \phi(x)$
- 動径基底関数回帰では、基底関数として「正規分布関数」を導入する。
- 基底関数「$\phi_h(x)$」は、以下のような図になる。
- これに適当な重み「$w=(w_{-H},w_{-H+1},\cdots,w_{H})^T$」を導入すると、関数「$y$」は以下のように表す事ができる。
$\Large \displaystyle \phi_h(x)=exp(-\frac{(x-\mu_h)^2}{\sigma~2})$
$\Large \displaystyle y=w^T \phi(x)=\sum_{h=-H}^{H}w_h exp(-\frac{(x-\mu_h)^2}{\sigma~2})$
- しかしながら、動径基底関数回帰には問題点があります。
- 例えば、グラフ上の「$-10\sim 10$」の間にメモリ1毎に動径基底関数を配置していくと、 以下のように多くの動径基底関数が必要になる「次元の呪い」が発生します。
- 「$x$」が2次元なら441個の動径基底関数が必要になる。 「$21^2=441$」
- 「$x$」が3次元なら9,261個の動径基底関数が必要になる。 「$21^3=9,261$」
- 「$x$」が5次元なら4,084,101個の動径基底関数が必要になる。「$21^5=4,084,101$」
- 動径基底関数回帰では、次元の呪いによって「$x$」の次元が大きくなるとすぐに、 現実的な計算可能リソースを超えてしまうのです。
4.2 ガウス過程
- 次元の呪いを解消する手段として、重みベクトル「$w$」について期待値をとって、 モデルから積分消去するというアイデアがある。
- 以下では、簡単のため「$\hat{y}=y$」と考える。
- 重みベクトル「$w$」が、平均が「$0$」で分散が「$\lambda^2 I$」のガウス分布から生成されたと考える。
- すると重みベクトル「$w$」の線形変換である「$y=\Phi w$」もガウス分布に従う。
$\Large \displaystyle \begin{pmatrix} \hat{y_1} \\ \hat{y_2} \\ \vdots \\ \hat{y_N} \end{pmatrix} = \begin{pmatrix} \phi_0(x_1) & \phi_1(x_1) & \cdots & \phi_H(x_1) \\ \phi_0(x_2) & \phi_1(x_2) & \cdots & \phi_H(x_2) \\ \vdots & \vdots & \ddots & \vdots \\ \phi_0(x_N) & \phi_1(x_N) & \cdots & \phi_H(x_N) \\ \end{pmatrix} \begin{pmatrix} w_0 \\ w_1 \\ \vdots \\ \vdots \\ w_H \end{pmatrix} $
$\Large \displaystyle \hat{y}=\Phi w$
$\Large \displaystyle w \sim \mathscr{N}(0,\lambda^2 I)$
- どんなN個の入力の集合「$(x_1,x_2,\cdots,x_N)$」についても、 対応する出力「$y=(y_1,y_2,\cdots,y_N)$」の同時分布「$p(y)$」が多変量ガウス分布に従う時、 「$x$」と「$y$」の関係はガウス過程(Gaussian process)に従う。
4.3 ガウス過程とカーネル
- ガウス過程回帰やサポートベクターなどで登場するカーネルとはデータポイント間がどれだけ似ているかを評価する指標である。
- 例えば下図の左のように丸、三角、四角の図形が並んでいたとする。
- この図形の中で似たもの同士をまとめていくつかのグループに分けてくれと言われたら、 あなたはどのようにグルーピングするか？
- 恐らく、多くの方は下図の右図の2パターンの分け方をするのではないだろうか？
- 上図のようにグルーピングする際に使った「色」や「形」といった指標こそがカーネルなのである。
- ガウス過程回帰を実行する際にもデータ間の類似度を評価するためにカーネルを使用するが、 どの指標で似ていると判断するかはカーネルによって異なる。
- 例えば、データ間のユークリッド距離の近さで判断するカーネルもあれば、 周期性や内積の大きさで判断するカーネルもある。
- カーネルとは何なのか？と言うことはなんとなく理解できたのではないでしょうか？
- ここからはカーネルの種類とその特徴について説明する。
- 今回は広く知られているカーネルとして以下の6つを紹介する。
- ガウスカーネル（RBF kernel）
- Matern カーネル（Matern kernel）
- 周期カーネル（Periodic kernel）
- 線形カーネル（Linear kernel）
- 多項式カーネル（Polynomial kernel）
- 指数カーネル（Exponential kernel）
- ガウスカーネルは上の数式で表されるカーネルである。
- 連続的な関数に使用され、無限次元の特徴空間に写像できる。
- データ間の類似性はガウス分布で評価し、滑らかな関数をモデル化する。
- ハイパーパラメータは、「$\sigma$」であり、データ間の類似性を評価する際に、どれだけ遠いデータを類似とみなすかを調整する。
- 「$\sigma$」が大きい場合、広い範囲のデータポイントを類似とみなす。
- 「$\sigma$」が小さい場合、近いデータポイントのみを類似とみなす。
- ガウシアンカーネルを使用するメリット、デメリットは以下の通りである。
- 〇メリット: 様々なデータパターンに適用可能であり、なめらかな関数をモデル化可能
- ☓デメリット: 適切なハイパーパラメータの選択が難しいことがあり、計算コストが高くなる可能性あり
$\Large \displaystyle K(x,x')=exp(-\frac{|x-x'|^2}{2\sigma^2})$
4.4 ガウス過程回帰モデル
4.5 ガウス過程のハイパーパラメータ推定
5.確率的生成モデル
5.1 自己相関関数とパワースペクトル
- ある波形と，それが一定値 だけずれたものがどれぐらい似ているかを表すものが，自己相関関数である．
- 時間関数で記述された信号の自己相関関数のフーリエ変換は，その信号のエネルギーあるいは電力の周波数軸上での分布を示す．
- 後でも述べるように時間関数では一意に記述できないような雑音等の不規則信号についても，自己相関関数は定義できる．
- 従って，確定信号と不規則信号を統一的に表現する上でも，自己相関関数とそのフーリエ変換は重要である。
- 自己相関関数「$\rho_f(t)$」は、次のように定義される。
$\Large \displaystyle \rho_f(t) =\int_{-\infty}^{\infty}f(t+\tau)f(\tau)d\tau$
- 二つの関数「$f_1(t)f_2(t)$」が与えらた時、以下の積分を重畳積分と呼ぶ。
- 時間域の重畳積分定理は、次のようになる
- 周波数域の重畳積分定理は、次のようになる
$\Large \displaystyle f(x)=\int_{-\infty}^{\infty}f_1(y)f_2(x-y)dy=f_1(x)*f_2(x)$
$\Large \displaystyle \int_{-\infty}^{\infty}f_1(\tau)f_2(t-\tau)d\tau \leftrightarrow F_1(\omega)*F_2(\omega) $
$\Large \displaystyle f_1(t)f_2(t) \leftrightarrow \frac{1}{2\pi}\int_{-\infty}^{\infty}F_1(y)F_2(\omega-y)dy$
- 周波数域の重畳積分定理から、Parsevalの等式を導く事ができる。
- 関数「$f(t)$」のフーリエ変換を「$F(\omega)$」とし、「$F(\omega)$」は次のように与えられると仮定する。
- 以下の「$E(\omega)$」を「エネルギースペクトル」と呼ぶ。
- 「$E(\omega)$」のフーリエ逆変換は、以下のようになる。
$\Large \displaystyle F(\omega) =A(\omega)e^{j\phi(\omega)}$
$\Large \displaystyle E(\omega) =A^2(\omega)$
$\Large \displaystyle E(\omega) =A^2(\omega)$
- 線形システムの入力「$f(t)$」に対する出力「$g(t)$」は、以下のようになる。
- 「$f(t)$」と「$g(t)$」の自己相関関数をそれぞれ「$\rho_f(t)$」と「$\rho_g(t)$」とすると、
- なる関係が成立する事を示そう。
- 「$f(t)$」と「$g(t)$」のフーリエ変換をそれぞれ「$F(\omega)$」と「$H(\omega)$」とすると、
$\Large \displaystyle g(t) =\int_{-\infty}^{\infty}f(\tau)h(t-\tau)d\tau=f(t)*g(t)$
$\Large \displaystyle \rho_g(t) =\rho_f(t)*h(t)*h(-t)$
$\Large \displaystyle g(t)\leftrightarrow F(\omega)*H(\omega)$
5.2 自己相関関数と周期性解析
以下の記事は、 「 自己相関関数と周期性解析」を引用しています。
- 製品開発においては不具合、不良をつぶしていかなければなりません。
- そういったものの中には設計に起因して規則性を持つ場合があります。
- 逆に言えば、規則性が理解できると設計との対応から不具合箇所を特定できる可能性があります。
- 今回は規則性を持つデータから周期性を見つけられる自己相関関数について述べます。
- 自己相関関数（ACF, auto correlated function）は時系列データの分析に用いられる手法です。
- 冒頭で述べたように製造工程の規則的な傾向や、設計に起因した規則性なども簡単に調べらます。
- 説明のために、まず「自己」という文字を取り除いた「相関関数」を考えてみます。
- 相関を確認するとき、よく用いられる指標は「相関係数」です。
- 相関関数というのは相関係数が「何か」の関数になっているということです。
- さらに自己相関であるからには自分自身との相関を意味します。
- 解析対象は一次元のデータです。
- そして時系列であったり、データアドレス順であったり、物理配置順であったりと、並びに意味があるデータです。
- このようなデータに対して自分自身との相関係数を考えます。
- 相関係数は(xi, yi)といったペアが対象がですから、データ並び順をずらした自分自身とのペアを考えます。
- このようにするとハッチングされたデータのペアの集合について相関係数が計算できます。
- 一つだけずらす必要もないので、ずらし量の分だけ沢山の相関係数が計算できることになります。
- 上図のずらし量をラグ（Lag）と呼びます。
- 自己相関関数は自己相関係数についてラグをパラメータとして関数としたものです。
- グラフ化するためにX軸にラグ、Y軸に自己相関関数をとったものを教科書的にはコレログラムと呼びます。
- Y軸は相関係数なので、コレログラムの縦軸は-1~+1になります。
- サンプルデータに自己相関関数を適用してみます。
- (1) 周期8+正規分布ノイズ小
- 並び順に意味を持った変数t=0~127に対し、x(t)としてtを8で割った余り(t mod 8)に標準偏差1の正規分布乱数を加えたものをテストデータとします。
- 上段のグラフにt - x(t)、下段に対象データのコレログラムを示します。
- 以下、同じ形式でいくつかサンプル解析を示します。
- 上段のグラフではなんとなく周期8の規則性が見える程度ですが、コレログラムの方でははっきりと8の倍数で正の相関が強いことがわかります。
- 周期8の依存性がある場合8の倍数周期も相関が高くなるので、周期的なピークの最大公約数が本当の周期性になります。
- (2) 周期8+正規分布ノイズ大
- 今度は(1)と同様で標準偏差2の正規分布乱数を加えたテストデータです。
- この例ではノイズが大きく上段のグラフでは周期性を読み取るのが辛くなってきますが、コレログラムは安定して8周期を捉えられています。
- (3) 周期8+周期16+正規分布ノイズ大
- 今度は(2)のデータにさらにtを16で割った余り÷2([t mod 16]/2)を加えてみます。
- 上段を見ると16周期に見えますが、コレログラムでは確かに8周期と16周期で正の相関が強そうです。
- 例(1)で倍数周期も見えてしまうと書きました。
- 通常ラグが大きくなると対となるデータが減少するので相関係数のバラツキも 大きくなります。
- しかし、この例では明らかに8周期より16周期の方が相関が高いです。 こういった例では16周期も意味があることが多いです。
- (4) まとめ
- 周期性を簡単に確認できる自己相関関数を紹介しました。
- この手法は信号解析や時系列解析で使われることが多いです。
-
しかし応用範囲は広く、製造工程や製品不良の周期のチェックなど周期的なものなら全て対象に出来ます。
（学問的にはウィーナー＝ヒンチンの定理に出会えたりして、ちょっとした感動を体験できます。） - またこちらでは順位相関係数を使って自己相関化数ライクなものをプロットしてみています。
- たとえノイズだらけのデータであっても、データ数が膨大な場合はコレログラムを図示するだけで色々気付く事があります。
- 盲目的に適用するだけでも、このアプローチを知っている人とそうでない人には大きな差がつくように思います。
6.ガウス過程の計算法
もくじ
6.1 カルマンフィルタリング問題
- 以下のようなスカラ時系列データ「$y(k)$」が、線形離散時間状態空間モデルで記述されると仮定する。
- ここで、「$x(k)$」は、n次元状態ベクトルである。
- 「$A$」は、n☓n行列である。
- 「$v(k)$」は、平均値「$0$」、分散「$\sigma_v^2$」の正規性白色雑音である。システム雑音とも呼ばれる。
- 「$w(k)$」は、平均値「$0$」、分散「$\sigma_w^2$」の正規性白色雑音である。観測雑音とも呼ばれる。
- 「$v(k)$」と「$w(k)$」は、お互いに独立であると仮定する。
- 状態空間モデルの係数行列・ベクトル「$(A,b,c)$」は、既知であると仮定する。
- この仮定は、時系列の状態空間モデルは、何等かの時系列モデリング法によって、既に得られているとする。
- また、係数行列・ベクトル「$(A,b,c)$」が時間によらず一定である時、上式で表される線形離散時間状態空間モデルは、 線形時不変（LTI: Linear Time-Invarient）システムとなり、その結果、得られる時系列「$y(k)$」は、定常確率過程になる。
- 以上の準備の元で、時刻kまでの時系列データ「$y(i),i=1,2,\cdots,k$」を用いて、 n次元状態変数ベクトル「$x(k)$」を推定する問題（フィルタリング問題）を考える。
- 特に、「$x(k)$」を何らかの意味で最もよく表す最適推定値「$\hat{x}(k)$」を最適フィルタを設計するものとする。
$\Large \displaystyle x(k+1)=Ax(k)+bv(k)$
$\Large \displaystyle y(k)=c^Tx(k)+w(k)$
- 最適制御と言う場合は、必ず評価関数（損失関数）を定義する必要がある。
- 評価関数（損失関数）の値が最小（又は最大）になる場合を最適制御と呼ぶ。
- 時系列データ「$y(i),i=1,2,\cdots,k$」に基づいて、状態「$x(k)$」のMSW（平均二乗誤差）の 最小値を与える推定値を見つける事を「カルマンフィルタリング問題」と言う。
- この推定値の事を「最小平均二乗誤差推定値（MMSE:Minimum Mean Square Error）」と呼ぶ。
6.2 逐次処理
6.3 時系列に対するカルマンフィルタ
- 時刻「$k$」における状態「$x(k)$」の推定値として、以下の二つの量を定義する。
- 事前推定値「$\hat{x}^{-}(k)$」
- 時刻「$k-1$」までに利用可能なデータに基づいた、時刻「$k$」における予測推定値を意味する。
- 事後推定値「$\hat{x}(k)$」
- 時刻「$k$」までに利用可能なデータに基づいた、即ち「$y(k)$」も用いた「$x」のフィルタリング時刻推定値。
- これが時刻「$k$」において求めるべき推定値である。
- カルマンフィルタでは、上図のように状態推定値は、時間更新されていく。
- 上図では、状態推定値（1次モーメント）「$\hat{x}(k-1)$」と共分散行列「$P(k-1)$」について示している。
- 図において、1時刻前の推定値「$\hat{x}(k-1)$」から現時点の事前推定値「$\hat{x}^{-}(k)$」までの処理を「予測ステップ」と呼ぶ。
- また、現時点の事前推定値「$\hat{x}^{-}(k)$」から現時点の事後推定値「$\hat{x}(k)$」までの処理を「フィルタリングステップ」と呼ぶ。
- 「最小二乗法」と同じように、ベイズ予測による以下のような事前推定値「$\hat{x}^{-}(k)$」と観測値「$y(k)$」に関する線形予測器を仮定する。
- 上式は、まさにベイズ推定であり、以下のような意味を持つ。
- 問題は「$G(k)$」と「$g(k)$」をどうやって決定するかである。
- ここでは、以下の事後状態推定誤差「$\tilde{x}(k)=x(k)-\hat{x}(k)$」に関する観測値「$y(k)$」との直交性を利用する。
$\Large \displaystyle \hat{x}(k)=G(k)\hat{x}^{-}(k)+g(k)y(k)$
事後推定値 = G(k)・事前推定値 + g(k)・観測値
$\Large \displaystyle \tilde{x}(k)=x(k)-\hat{x}(k)$
$\Large \displaystyle E[\tilde{x}(k)y(i)]=0,\ \ i=1,2,\cdots,k-1$
6.4 数値シュミレーション例
6.5 システム制御のためのカルマンフィルタ
以上
7.ガウス過程の適用
もくじ
7.1 LTIシステムの表現
- システムへの入力信号を「$u(k)$」、出力信号を「$y(k)$」とすると、
- なる「畳み込み積分」が成り立つ。ただし、「$g(k)$」は、システムのインパルス応答である。
$\Large \displaystyle y(k)=\sum_{l=0}^{\infty}g(l)u(k-l)$
- システムへの入力信号を「$u(t)$」、出力信号を「$y(t)$」とすると、
- なる「畳み込み積分」が成り立つ。ただし、「$g(t)$」は、システムのインパルス応答である。
$\Large \displaystyle y_t=\int_0^t g(\tau)u(t-\tau)d\tau=g(t)*u(t)$
7.2 入出力関係
- 上記の離散時間LTIシステムをz変換すると、以下の式が得られる。
- 但し、「$u(z)$」「$y(z)$」は、それぞれ「$u(k)$」「$y(k)$」のz変換であり、以下のように定義される。
- 但し、「$G(z)$」は、「$g(k)$」のz変換であり、以下のように定義され、離散時間LTIシステムの伝達関数と呼ばれる。
- 一般に、離散時間LTIシステムは、以下のような線形定係数差分方程式で記述できる。
- 上式を初期値をゼロとして、z変換すると、以下のようになる。
- 従って、以下のような伝達関数「$G(z)$」が得られる。
$\Large \displaystyle y(z)=G(z)u(z)$
$\Large \displaystyle u(z)=\sum_{l=0}^{\infty}u(k)z^{-k}$
$\Large \displaystyle y(z)=\sum_{l=0}^{\infty}y(k)z^{-k}$
$\Large \displaystyle G(z)=\sum_{l=0}^{\infty}g(k)z^{-k}$
$\Large \displaystyle y(k)+a_1 y(k-1)+\cdots+a_n y(k-n)=b_1 u(k-1)+b_2 u(k-2)+\cdots+b_mu(k-m)$
$\Large \displaystyle (1+a_1 z^{-1}+\cdots+a_n z^{-n})y(z)=(b_1 z^{-1}+b_2 z^{-2}+\cdots+b_m z^{-m})u(z)$
$\Large \displaystyle G(z)=\frac{y(z)}{u(z)}=\frac{b_1 z^{-1}+b_2 z^{-2}+\cdots+b_m z^{-m}}{1+a_1 z^{-1}+\cdots+a_n z^{-n}}$
7.3 ARXモデル
- シフトオペレータ「$q$」は、以下のような性質を持つ。
$\Large \displaystyle qy(k)=y(k+1)$
$\Large \displaystyle q^{-1}y(k)=y(k-1)$
- ARXモデルは、ARモデルに外生（因）入力(eXogenous)が加わったモデルである。
- 従って、ARXモデルの式は、以下のように定義される。
- ここで、「$y$」は出力、「$u$」は入力、「$w$」は誤差やノイズを意味し、「$k$」は離散時間系におけるサンプル点を意味する。
- 「$q$」はシフトオペレータと呼ばれるもので、「$z$」変換における演算子「$z$」と意味は同じである。
- ARXモデルのブロック図は下記のようになる。
- すなわちARXモデルは、コントローラによって決定される操作量「$u$」とプラントからの出力「$\frac{w(k)}{A(z)}$」を利用したシステム同定であり、ノイズは正確な同定を阻害する因子と言える。
- 非常に分かり易いですね。
$\Large \displaystyle A(q)y(k)=B(q)u(k)+w(k)$
以上
8.ガウス過程による教師なし学習
もくじ
8.1 カルマンフィルタ
- 測定には誤差がつきものである。誤差を減らす策としては、同一条件で独立な測定を何度も行う事である。
- 「独立同一分布（i.i.d.）」と言う前提の元でn回測定して平均すれば、大数の法則により標準偏差が「$\frac{1}{\sqrt{2}}$」になる。
- しかし「止まった的を射る」には、この方法でいいのだが、「動いている的を射る」には、この方法では上手く行かない。
- 時刻kにおける的の位置を、実数値の確率変数「$x(k)$」で表す事にする。
- 我々は、的の位置である「$x(k)$」を当てたいが、直接知る事は出来ないとする。
- 従って、代わりに観測誤差「$w(k)$」を含む観測値「$y(k)$」を得る事が出来るとする。
- 次に、現在の的の位置である「$x(k)$」」から、次の時間の的の位置である「$x(k+1)$」を予想する。
- 「k」と「k+1」が短い時間差であれば、「$x(k+1)$」は「$x(k)$」」から少しだけの差「$v(k)$」を足した位置と予想できる。
- 従って、以下のような式の確率モデルを作る事が出来る。
- ただし、以下のような確率分布（偏差は定数）を仮定し、お互いに独立であるとする。
- つまり、カルマンフィルタの問題設定は、観測値「$y(1),y(2),\cdots, y(k)$」からし、状態変数(位置$x(k)$)を推定する事になる。
$\Large \displaystyle x(k+1)=x(k)+v(k)$
$\Large \displaystyle y(k)=x(k)+w(k)$
$\Large \displaystyle x(0)\sim N(0,\sigma_0^2)$
$\Large \displaystyle v(k)\sim N(0,\sigma_v^2)$
$\Large \displaystyle w(k)\sim N(0,\sigma_w^2)$
- 信号「$x$」は、以下の観測式で観測されるものとする。
- ただし、「$y$」は、観測値であり、「$w$」は、正規分布「$N(0,\sigma^2)$」に従う観測雑音とする。
- 推定問題は、雑音に汚された観測値「$y$」から、信号「$x$」の推定値「$\hat{x}$」を求める事である。
- さて、線形推定則を仮定すると、信号「$x$」の推定値「$\hat{x}$」は、以下のような観測値「$y$」の1次式で表される。
- 以下のような推定誤差(error estimate)を定義し、この平均二乗誤差（MSE: Mean Square Error）が 最小になるように、線形推定則のパラメータ「$\alpha$」と「$\beta$」を決定する問題を考える。
- 以下では、「①推定誤差の平均がゼロにする」と「②推定誤差の分散を最小にする」と言う 二つの条件を使って、未知数「$\alpha,\beta$」を決定する。
$\Large \displaystyle y=cx +w$
$\Large \displaystyle \hat{x}=f(y)=\alpha y + \beta$
$\Large \displaystyle e=x-\hat{x}$
- この条件を計算すると、
- となるから、
- の時、推定値はバイアスを持たない。この時「$\hat{x}$」は不偏推定値と呼ばれる。
$\Large \displaystyle E[e]=E[x-\hat{x}]=E[x-\alpha y - \beta]$
$\Large \displaystyle =E[x-\alpha (cx+w) - \beta]$
$\Large \displaystyle =(1-\alpha c)\bar{x}-\alpha\bar{w}-\beta=0$
$\Large \displaystyle \hat{\beta}=(1-\alpha c)\bar{x}-\alpha\bar{w}=\bar{x}-\alpha(c\bar{x}+\bar{w})$
- 推定値の平均二乗誤差（MSE）は、以下のようになる。
- 上式を「$\alpha$」に関する2次式に書き直して平方完成すると以下のようになる。
- いま、
- とおくと、上式は、以下のように変形できる。
- これより、
- の時、推定誤差分散は、最小値、
- をとる事は分かる。
- この時、「$\hat{x}$」は、「最小分散推定値」と呼ばれる。
$\Large \displaystyle E[\{e-E(e)\}^2]=E[[(x-\alpha y-\beta)-\{(1-\alpha c)\bar{x}-\alpha \bar{w}-\beta \}]^2]$
$\Large \displaystyle =E[\{ x-\alpha(cx+w)-(1-\alpha c)\bar{x}+\alpha \bar{w}\}^2]$
$\Large \displaystyle =E[\{ (1-\alpha c)(x-\bar{x})-\alpha(w-\bar{w}) \}^2]$
$\Large \displaystyle =E[(1-\alpha c)^2(x-\bar{x})^2+\alpha^2(w-\bar{w})^2 -2\alpha(1-\alpha c)(x-\bar{x})(w-\bar{w}) ]$
$\Large \displaystyle =(1-\alpha c)^2E[(x-\bar{x})^2]+\alpha^2E[(w-\bar{w})^2] -2\alpha(1-\alpha c)E[(x-\bar{x})(w-\bar{w})] $
$\Large \displaystyle =(1-\alpha c)^2\sigma_x^2+\alpha^2\sigma_w^2$
$\Large \displaystyle E[\{e-E(e)\}^2]=(1-2\alpha c+\alpha^2 c^2)\sigma_x^2+\alpha^2\sigma_w^2$
$\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)\alpha^2 - 2c\sigma_x^2 \alpha + \sigma_x^2$
$\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)[\alpha^2 -\frac{2c\sigma_x^2}{c^2 \sigma_x^2 +\sigma_w^2}\alpha+\frac{\sigma_w^2}{c^2 \sigma_x^2 +\sigma_w^2} ]$
$\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2})^2 - \frac{c^2 \sigma_x^4}{c^2 \sigma_x^2 + \sigma_w^2} + \sigma_x^2$
$\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_x^2}{c^2 \sigma_x^2 + \sigma_w^2})^2 - \frac{\sigma_x^2 \sigma_w^2}{c^2 \sigma_x^2 + \sigma_w^2}$
$\Large \displaystyle =(c^2 \sigma_x^2 +\sigma_w^2)(\alpha -\frac{c\sigma_w^{-2}}{c^2 \sigma_w^{-2} + \sigma_x^{-2}})^2 - \frac{1}{c^2 \sigma_w^{-2} + \sigma_x^{-2}}$
$\Large \displaystyle \sigma^2=\frac{1}{c^2 \sigma_w^{-2} + \sigma_x^{-2}}$
$\Large \displaystyle E[\{e-E(e)\}^2]=(c^2 \sigma_x^2+\sigma_w^2)(\alpha - c\sigma_w^{-2}\sigma^2)^2+\sigma^2$
$\Large \displaystyle \alpha=c\frac{\sigma^2}{\sigma_w^2}$
$\Large \displaystyle E[\{e-E(e)\}^2]=\sigma^2$
- 以下のようなスカラ時系列データ「$y(k)$」が、線形離散時間状態空間モデルで記述されると仮定する。
- ここで、「$x(k)$」は、n次元状態ベクトルである。
- 「$A$」は、n☓n行列である。
- 「$v(k)$」は、平均値「$0$」、分散「$\sigma_v^2$」の正規性白色雑音である。システム雑音とも呼ばれる。
- 「$w(k)$」は、平均値「$0$」、分散「$\sigma_w^2$」の正規性白色雑音である。観測雑音とも呼ばれる。
- 「$v(k)$」と「$w(k)$」は、お互いに独立であると仮定する。
- 状態空間モデルの係数行列・ベクトル「$(A,b,c)$」は、既知であると仮定する。
- この仮定は、時系列の状態空間モデルは、何等かの時系列モデリング法によって、既に得られているとする。
$\Large \displaystyle x(k+1)=Ax(k)+bv(k)$
$\Large \displaystyle y(k)=c^Tx(k)+w(k)$
- まず、k=1の場合を考えて見る。「$y(1)$」を見て、「$x(1)$」を推定するために、両者の同時分布を求める必要がある。
- 設定より、以下の式が成立する。
- 行列を使うと、以下のように書ける。
- 右辺のベクトル「$(w(1),v(0),x(0))^T$」は、前提から三次元正規分布に従う。
- 従って、左辺のベクトル「$(y(1),x(1))^T$」も、二次元正規分布に従う。
$\Large \displaystyle x(1)=x(0)+v(0)$
$\Large \displaystyle y(1)=x(1)+w(1)=x(0)+v(0)+w(1)$
$\Large \displaystyle \begin{pmatrix} y(1) \\ x(1) \end{pmatrix} = \begin{pmatrix} 1 & 1 & 1 \\ 0 & 1 & 1 \end{pmatrix} \begin{pmatrix} w(1) \\ v(0) \\ x(0) \end{pmatrix} $
- 線形離散時間定常確率システムは，下図のように表現される。
- 線形離散時間定常確率システムの状態方程式は，次のようになる。
- 「$x(t+1)$」は、1期先の状態変数である。「$x(t)$」は、現在の状態変数である。
- 「$v(k)$」は、「$v(k)\sim N(0,\sigma_v^2)$」のシステム雑音（神の見えざる手）である。
- 「$A\in R^{n\times n}$」および「$b\in R^{n}$」である。
- 線形離散時間定常確率システムの観測方程式は，次のようになる。
- 「$x(t)$」は、現在の状態変数である。
- 「$w(k)$」は、「$w(k)\sim N(0,\sigma_w^2)$」の観測雑音である。
- 「$C\in R^{n}$」である。
- カルマンフィルタでは以下の補正の式を使います。
- 補正後の状態 ＝ 補正前の状態 ＋ カルマンゲイン ×（本物の観測値－予測された観測値）
- カルマンゲイン ＝ 状態の予測誤差の分散 ÷（状態の予測誤差の分散＋観測方程式のノイズの分散）
$\Large \displaystyle x(k+1)=Ax(k)+bv(k)$
$\Large \displaystyle y(k)=c^Tx(k)+w(k)$