Variational AutoEncoder( VAE )

近年、ディープラーニング業界、はたまた画像処理業界ではGANとVAEの2つの技術で話題が持ちきりとなりました。GANについては前回解説しました。今回はそんなVAEのアルゴリズムについて解説をしていきます。GAN同様に最後にPyTorchによる実装例を紹介していきます。

VAEとは

VAEは情報を圧縮して圧縮した情報から元の情報に戻す、というような仕組みをもった、AE(Auto encoder)と言われるものの一種です。AEはただ単にデータの圧縮と再構築をするだけでしたが、VAEは確率モデルの理論を導入しています。VAEは確率分布を使ったモデルということは、未知のデータを確率的に作成できることになります。VAEはGenerative model(生成モデル)と言われています。トレーニングデータからProbability Density Function (PDF)を推定するモデルであるためす。 

GANと違って何ができるのか

VAEの一つの特徴は次元削減です。次元削減といえばPCASVDなどを思い浮かべるかもしれません。VAEは同様にLatent space(潜在空間)に次元圧縮し、またそこから復元するということをしています。

PCAのように圧縮した次元において直行している必要はありません。この辺も普通の次元圧縮とのアプローチとは異なっています。

赤はAEの軸。青はPCAの軸。PCAは軸同士は直行しているが、AEに関しては必ずしも直行している必要はない。

AEとVAEは違いがあまり無く感じますが、前述の通り確率を用いたところに違いがあります。AEは点(Single point)としてデータを潜在空間にマッピングしていましたが、VAEでは潜在空間に落とし込むときにガウス分布に従って落とし込ませます。デコーダーはその点を拾ってデコードするため、見方を変えれば、確率的にサンプリングしてzを拾っている、ということになるのです。このzをデコードして元の入力に近づけるようにするのです。

さて、VAEは確率モデリングと言われます。確率モデリングとは、例えば何かのデータxが分布していたとして、その分布を確率で表現するというモデルです。あるデータxが正規分布に従っていたとしたら、データxの確率分布は正規分布となり、P(X)=正規分布( μ, σ )と言う感じで表現します。P(X)と書くと思わずPは関数と勘違いするかもしれませんが、確率を表すものなので注意が必要です。

数式

VAEはデータを確率モデル化をすることを目標とします。データをX、その確率分布をP(X)と定義すれば、P(X)を見つけること、すなわちP(X)の最大化がVAEの目的です。そしてP(X)に関しては次のように表現することができます。[1]

P(X) = \int P(X \vert z) P(z) dz = \int P(X,z)

上記の式はXが入力画像と考え、その画像を表現する潜在ベクトルzを少なくとも1つ見つけたいという意味があります。

また、事前分布P(z)から\{z_i\}_{i=1}^nをサンプルした際に、P(X)を次のように近似できます。

P(X) \approx \frac{1}{n}\sum_{i=1}^n P(x|z_i)

P(X)を求めればいいのですが、一筋縄には行きません。それはXが高次元であることもそうですが、その確率を求めるのにはすべてのzの空間を舐め回すような非常にたくさんのサンプリングが必要なだけでなく、組み合わせをチェックするような処理が必要となり、総当りでP(X)を求めることは現実的ではありません。

そのため、P(X)を求めるには別のアプローチを考えます。もし事後確率のP(z|X)がもとめられれば、p(x|X)=\int P(x|z)p(z|X) dz より未知の画像xを作り出すような確率分布を求めることができそうです。そうなると 事後確率のP(z|X)を求めるだけとなり、簡単に思えますが事後確率はベイズの定理により次の式になります。

P(z|X) = P(X|z)\cfrac{P(z)}{ P(X)}

よく見るとあの厄介なP(X)が分母にあります。やはり求めることができません。

ここで、諦めずに別のアプローチで P(z|X) を求めていくようにしていきます。

P(z|X) を求める

P(z|X)を求めるにあたり少し工夫します。ここでVAEの名前の由来ともなっていますが、 Variational Inference(VI) という手法を使います。P(z|X)を推定するためにQ(z|X)という分布を考えます。Q(z|X)はP(z|X)の近似で、ガウス分布などの簡単な分布関数の組み合わせによりP(z|X)を近似します。

関数の組み合わせにより、近似していく。青の点線へ緑の分布関数で近似していく様子。 https://towardsdatascience.com/bayesian-inference-problem-mcmc-and-variational-inference-25a8aa9bce29

どれだけ似ているのか、という指標については類似尺度としてKLダイバージェンスを使います。同じ分布であれば0に、異なれば値が大きくなっていくような関数です。KLダイバージェンス次のように書くことが出来ます。

KL(P||Q) = E_{x \sim P(x)} log \cfrac{P(x)}{Q(x)}=\int_{-\infty}^{\infty}P(x)\ln \cfrac{P(x)}{Q(x)}dx

ここから式が多く出てくるので、見やすくするために単純にP(z|X )をP, Q( z|X)をQと表示します。Qの近似は次のようにかくことができます。

\begin{aligned} D_{KL}[Q\Vert P] &= \sum_z Q \log (\cfrac{Q}{P}) \\ &=E [ \log (\cfrac{Q}{P}) ] \\ &= E[\log Q - \log P] \end{aligned}

最後はlogの変換公式により、割り算を引き算にしただけです。さてP である P(z|X ) は P(z|X) = P(X|z)P(z) / P(X) のように表現できたので、先の式に当てはめてみましょう。

\begin{aligned} D_{KL}[Q\Vert P] &= E[\log Q - (\log (P(X \vert z) \cfrac{P(z)}{P(X)})] \\ &= E[\log Q- \log P(X \vert z) - \log P(z) + \log P(X)] \end{aligned}

ここでzに関する期待値でくくっている項に注目するとP(X)があります。P(X)はzに関係ありません。そのため括弧の外に出すことができます。

\begin{aligned} D_{KL}[Q\Vert P] &= E[\log Q- \log P(X \vert z) - \log P(z)] + \log P(X) \\ D_{KL}[Q\Vert P] - \log P(X) &= E[\log Q- \log P(X \vert z) - \log P(z)] \end{aligned}

ここからがトリッキーなのですが、右辺に注目するともう一つのKLダイバージェンスを見つけることができます。先の式を両辺にマイナス倍します。

\begin{aligned} \log P(X) - D_{KL}[Q\Vert P] &= E[-\log Q + \log P(X \vert z) + \log P(z)] \\ &= E[\log P(X \vert z) -( \log Q - \log P(z))] \\ &= E[\log P(X \vert z)] -E[( \log Q - \log P(z))] \\ &= E[\log P(X \vert z)]- D_{KL}[Q\Vert P(z) ] \end{aligned}

なんとPとQを近づけようとした結果、本来の目的であったP(X)に関する式がでてきてしまいました。PとQを縮約して記載していたので、正しく展開してかいてみます。

\log P(X) - D_{KL}[Q(z \vert X) \Vert P(z \vert X)] = E[\log P(X \vert z)] - D_{KL}[Q(z \vert X) \Vert P(z)]

本当の目的はP(X)の最大化でした。結局ですが、これを最終的なVAEの目的の関数として定義することになります。

最終的に得た式は非常に興味深い構造になっています。

  1. Q(z|X) はデータXを受取り、潜在空間にzを投影
  2. zは潜在変数
  3. P(X|z)は潜在変数zからデータ生成

つまりQ(z|X)はエンコーダ、zは潜在変数(エンコードされたデータ),そしてP(X|z)はデコーダです。まさにオートエンコーダーそのものとなりました。

算出された式を観察して、左辺と右辺がありますが、左辺にあるlog(p)最大化することが目的でした。そのためには右辺を最大化していくこと、つまり E[\log P(X \vert z)] を大きくして D_{KL}[Q(z \vert X) \Vert P(z)]を小さくしていくことで、左辺は大きくなっていきます。

そのため、今度は目的を右辺の最大化に絞っていきます。

右辺を最大化する。

右辺には以下の2つの式があります。

  1. E[\log P(X \vert z)]
  2. D_{KL}[Q(z \vert X) \Vert P(z)]

右辺=数式1-数式2、ですので数式1を最大化、数式2は最小化していけば、右辺は大きくなり目的が達成されます。

まず 数式 1についてですが、よく見ると潜在変数zを受取りXを生成する教師つきの学習そのものです。ですので、学習によりなんとかなりそうです。

さて、厄介なのが 数式 2です。ここで一つの仮定を起きます。P(z)は正規分布N(0, 1)と仮定するのです。そして、Xからzを生成する分布もパラメータ\mu(X), \sigma(X)付きの正規分布となります。 平均と分散はXを中心としたという意味です。そして、KLダイバージェンスは次のように表されます。

D_{KL}[N(\mu(X), \Sigma(X)) \Vert N(0, 1)] = \frac{1}{2} \, \left( \textrm{tr}(\Sigma(X)) + \mu(X)^T\mu(X) - k - \log \, \det(\Sigma(X)) \right)

kはガウシアンの次元数、traceは対角要素の和を表します。そして、detは対角要素の積\det \left({\mathbf A}\right) = \prod_{i \mathop = 1}^n a_{ii}です。

導出に関しては、ここでは重要でないため割愛しますが最終的には次のようになります。(導出に関して興味ある方は[2]を参照してください。)

D_{KL}[N(\mu(X), \Sigma(X)) \Vert N(0, 1)] = \frac{1}{2} \sum_k \left( \exp(\Sigma(X)) + \mu^2(X) - 1 - \Sigma(X) \right)

この項目を実装するときにはロス関数の中で利用します。実際との差分を計算するためですね。

全ての要素の説明ができました。あとはやるだけです。

実装

エンコーダーとデコーダーの実装は次のとおりです。

エンコーダー

class Encoder( nn.Module ):
    def __init__( self ):
        super().__init__()
	self.common = nn.Sequential(
            nn.Linear( 784, 400 ),
            nn.ReLU(),
            )
	self.model1 = nn.Sequential(
            self.common,
            nn.Linear( 400, 20 )
            )
        self.model2 = nn.Sequential(
            self.common,
            nn.Linear( 400, 20 )
            )
    def forward( self, img ):
	img_flat = img.view( img.size( 0 ), -1 )
        return self.model1( img_flat ), self.model2( img_flat )

デコーダー

class Decoder( nn.Module ):
    def __init__( self ):
	super().__init__()
	self.model = nn.Sequential(
            nn.Linear( 20, 400 ),
            nn.ReLU(),
            nn.Linear( 400, 784 ),
            nn.Sigmoid(),
            )
    def forward( self, z ):
        return self.model( z )

エンコーダーとデコーダーを用いいたVAEを次のように実装していきます。

VAE

class VAE( nn.Module ):
    def __init__( self ):
        super().__init__()
        self.encoder = Encoder()
	self.decoder = Decoder()

    def _reparameterization_trick( self, mu, logvar ):
        std = torch.exp( 0.5 * logvar )
        eps = torch.randn_like( std )
        return mu + eps * std

    def forward( self, _input ):
        mu, sigma = self.encoder( _input )
	z         = self._reparameterization_trick( mu, sigma )
        return self.decoder( z ), mu, sigma

説明ではzの取得は正規分布でのサンプリングを仮定しました。ところが実際にサンプリングするとバックプロパゲーションが出来ないため学習が出来ません。そこでreparameterization trickというのを用いいます。zを次の式で近似します。

z = \mu(X) + \Sigma^{\frac{1}{2}}(X) \, \epsilon

where, \epsilon \sim N(0, 1)

左図はZを表現するために本来のサンプリングで実装した図。左ではバックプロパゲーションができない。そのため、足し算と掛け算に依る表現に正規分布に従うノイズの足しこみを行い誤差伝搬を可能にする。[3]

最後に損失関数とVAEを使ったコードは次のとおりです。

# Kingma and Welling. Auto-Encoding Variational Bayes. ICLR, 2014                                                                                                                                                                 
# 入力画像をどのくらい正確に復元できたか?                                                                                                                                                                                        
def VAE_LOSS( recon_x, x, mu, logvar ):
    # 数式では対数尤度の最大化だが交差エントロピーlossの最小化と等価                                                                                                                                                              
    BCE = F.binary_cross_entropy(recon_x, x.view(-1, 784), size_average=False)
    # 潜在空間zに対する正則化項. # P(z|x) が N(0, I)に近くなる(KL-distanceが小さくなる)ようにする                                                                                                                               
    KLD = -0.5 * torch.sum(1 + logvar - mu.pow(2) - logvar.exp())
    return BCE + KLD

def main():
    epoch_size = 50
    vae = VAE()
    vae.cuda()
    Tensor = torch.cuda.FloatTensor
    dataloader=get_dataloader()
    optimizer = torch.optim.Adam( vae.parameters(), lr=1e-3 )

    for epoch in range( epoch_size ):
        for i, ( imgs, _ ) in enumerate(dataloader):
            optimizer.zero_grad()
            real_images          = Variable( imgs.type( Tensor ) )
            gen_imgs, mu, logvar = vae( real_images )
            loss                 = VAE_LOSS( gen_imgs, real_images, mu, logvar ).cuda()
            loss.backward()
            optimizer.step()

VAE_LOSSでは得た画像が目的とした式と同じになるような差分を計算しています。

長くなりましたがこれがVAEの全貌です。今回のソースコードはGithubにあげてあります。

https://github.com/octopt/techblog/blob/master/vae/main.py

参考文献

  1. https://en.wikipedia.org/wiki/Law_of_total_probability
  2. https://wiseodd.github.io/techblog/2016/12/10/variational-autoencoder/
  3. https://towardsdatascience.com/understanding-variational-autoencoders-vaes-f70510919f73

おすすめの記事

最小二乗法による重回帰分析

重回帰分析の値を最小二乗法でえる方法を示す。早速だが、次の式に適用することで目的の係数群が手に入る。

\mathbf{w} = (\mathbf{X}^{T}\mathbf{X})^{-1} \mathbf{X}^{T}\mathbf{y}

上記x, y, wは次の通りの関係を持っている。

y_i=a_{0i} + x_{0i}w_{0i} + x_{1i} w_{1i} + ... x_{ni} w_{ni}+ f

aは切片、fは残差と言われる。ここでyはm個あると仮定する。回帰分析のところで説明したとおり、切片aと重みwベクトルを求めるのが回帰分析の目的となる。最初の\mathbf{X}については次のようになる。

\mathbf{X} = \begin{bmatrix} \mathbf{x}_0 & \mathbf{x}_1 & … &\mathbf{x}_n \end{bmatrix}

\mathbf{x}_0 = \begin{bmatrix} x_{00} \\ x_{01} \\ ... \\ x_{0m} \end{bmatrix}

となっている。

上記ではやや理解がしづらいと思うので、実際に具体例をしめす。なお、他のサイトでは切片を求めない場合もある。当サイトでは、今回は切片を考慮に入れて説明しよう。

LSMを用いた重回帰分析の仕方(計算方法具体例)

具体的に説明のためにあるサイト[2]から取得したデータを用いる

営業所売上額広告費販売員数
A85006
B95008
C1370010
D1140013
E1480011
F17120013
G?130014

上記で営業所Gの売上高を予想したいというのが今回の目的である。(そのために重みベクトルと切片を求める。)上記の表で例えば営業所Aの売上高は 8 = a_0 + 500 * w_0 + 6 * w_1 ということになり、最小二乗法を使うことで未知の切片のa_0と重みである\mathbf{w}が求まる。

まずは上記表を元にa0, x0, x1のベクトル、並びにyのベクトルを定義する。

import numpy
def main():
    y = numpy.array( [[8, 9, 13, 11, 14, 17]]).T
    a0 = numpy.array( [[ 1, 1, 1, 1, 1, 1]] ).T # 切片が不要な場合はここは不要
    x1 = numpy.array( [[ 500, 500, 700, 400, 800, 1200]] ).T
    x2 = numpy.array( [[ 6, 8, 10, 13, 11, 13]] ).T

a0に関しては全て1にする。切片であるために1以外の数は適切ではない。そして\mathbf{X}は次のように定義される。

X = numpy.hstack( [a0, x1, x2] ) # 切片が不要であれば最初のa0をなくす。
# [[   1  500    6]
# [   1  500    8]
# [   1  700   10]
# [   1  400   13]
# [   1  800   11]
# [   1 1200   13]]

肝心のXは上記のように出力される。なお、切片が不要な場合は最初の「1」がなくなった行列となる。さて、定義できたのであとは最初に示した式通りに書いていく。

XTX = X.T.dot(X) # X^T.X                                                                                                                                                                                                      
invXTX = numpy.linalg.inv( XTX )
result = ( invXTX.dot( X.T ) ).dot( y )
print( result )
# [[1.14814815]
#  [0.00786008]
#  [0.53909465]]

出力値として、a_0=1.14814815, x_0=0.00786008, x_1=0.53909465となった。

この値を元にGの営業所の売上予想は 1.14814815 + 1300 \times 0.00786008 + 14 \times 0.53909465=18.9135772となる。

参考文献

  1. https://k-san.link/linear-regression-01/
  2. https://istat.co.jp/ta_commentary/multiple
  3. https://datachemeng.com/wp-content/uploads/ordinaryleastsquares.pdf

おすすめのサイト

数量化理論Ⅳ類(MDS)

数量化理論Ⅳ類はMDS,多次元尺度構成法とも呼ばれます。対象物の類似度を測る手法です。

あるアイテムが似ている、似ていないというのを平面上にプロットすることができます。平面上だけでなくても、直線、或いは3次元上でプロットできますが、多くは視覚的にわかりやすい2次元にプロットする例が多いです。なお、「似ている程度」「似 ていない程度」の量のことを数量化理論Ⅳ類では親近性と呼びます。

MDSの説明ですが、以下のドキュメントが非常にわかりやすく書かれています。導出が気になった方は、是非次のドキュメントに目を通してください。

http://www.page.sannet.ne.jp/yo-skmt/MultiEx/No13-15MU-Text02.pdf

女性ファッション雑誌をMDSで可視化

MDSの利用シーンとしての具体例を考えます。例えば女性専門雑誌がいくつかあるとします(VIVI, an anなど)。女性向けファッション雑誌の内容は男性からするとどれも一緒に見えますが、女性から見ればファッションの種類によって雑誌がかわります。では具体的にどの雑誌とどの雑誌が似ているのか、という事を男性が知りたいときにはどうすればいいのでしょうか。(男性は例え本の中身を見てもさっぱりですよね?)

まさにこんなときがMDSの出番と言えるでしょう。MDSを使えば次のようにマッピングできます。

上記画像は https://u-site.jp/research/methods/mds/ より。

雑誌の中身を見て、これとこれが似ているなど判断しなくてもMDSによりわかりやすく可視化することができます。

その他にも、例えばM1グランプリの採点から、採点者の傾向や、お笑い芸人のタイプなどを可視化出来ます。採点表はネットに公開されているので、興味があったらMDSで実験してみたら面白いかもしれません。

MDSの計算方法

対象物間に「距離」という概念を持ち込む必要があります。対象物間どうしの「距離」の概念を何らかの方法で定義できれば、 固有値分解を実施することでMDSが完了します。距離とは何か、ということを疑問に思った方は距離の公理を参照ください。大抵の場合は馴染み深い距離(ユークリッド距離)を利用するケースが多いです。

ちなみに、どのような行列に対しての固有値分解(あるいは特異値分解)をすることで分析手法の名前が変わってきます[2]

  1. 相関行列の固有値分解⇒主成分分析or因子分析 (数量化理論III類)
  2. 分散比行列の固有値分解 ⇒ 判別分析 (数量化理論II類)
  3. 頻度行列の特異値分解 ⇒コレスポンデンス分析(数量化理論III類)
  4. 距離行列の固有値分解 ⇒ MDS (数量化理論IV類)

MDSは対象間同士の距離行列さえ求めてしまえば、あとはやる(主成分分解する)だけで可視化できるというのはある意味興味深いです。

距離行列の作成

実際の例で計算をして見ます。今、A,B,Cの三人がいたとして、それぞれの互いにどう思っているかを評価してもらいます。(親近性の調査)

AはBはC は
AのことNaN33
B のこと 4 NaN 3
C のこと 32 NaN

上記のように、互いに評価をしています。なおこの行列を親近性行列と呼びます。

MDSとは、例えば3人を親密な人は距離が近く、そうでない人同士が遠くなるように位置を求めることです。2次元でも、3次元でも考えられますが、計算できる次元数は{対象者−1}次元となります。なお、親近性行列の要素は e_{11}, e_{12}, e_{13}, e_{21} ... e_{32}, e_{33}のように、 e_{ij}として表すことにします。

1次元の直線でA,B,Cの親近性を表現すること考えてみます。A, B, Cのそれぞれの位置をx1, x2, x3とします。ここで、位置x1, x2, x3の制約条件を次のとおりに定めます。

  • x1, x2, x3の平均は0
  • x1, x2, x3の分散は1

そして関数Qというものを定義します。関数Qは次のような定義です

m=特定の2人の親近性\times 特定の 2人の距離の2乗

つまり次のようになります。

m=e_{ij}\times (x_{i} - x_{j})^2

全員分のmの総和にマイナスしたものをQとして定義します。

Q=-\sum{ m_k }

Qの最大化を行うと、親近性の大きい個体の距離は短く、親近性の小さい個体の距離がなくなるようになります。なお、距離の定義は今回は線分で見るのでユークリッド距離を使いますが対象により変わります。要は距離の公理の「距離」が定義できるものであれば何でも良いです。

さてこれらを用いて計算ですが、細かい導出が坂元保秀教授により大変素晴らしくまとめられています。[1]

数式が多くなるので、結論を書きますと親近性行列から行列Hを定義します。行列Hとはh_{12} = e_{12} + e_{21}という形で、対角線同士を足し合わせる行列となります。対角線の要素[ latex ] e_{11}, e_{22}, e_{33} [/latex]については、e_{11} = -(e_{12}+e_{13}) \times 2 のように、自分意外の行を全て足しこみ、二倍してマイナスにした数字です。以上の操作を図にすると次のとおりです。

行列Hの作り方。親近性行列と、その天地を足す。加えて、親近性行列の行をすべて足しこんだものを2倍してマイナスをかけたものを対角の要素とする

行列Hがでたあとは、固有値分解をすることでx1,x2,x3が算出されます。いま、Pythonを用いて行列Hを固有値分解してみましょう。

import numpy as np
import math

a = [-12, 7, 6]
b = [7, -14, 5]
c = [ 6, 5, -5]
mat = np.array( [ a, b, c ] )
e_value, e_mat = np.linalg.eig( mat )
x1, x2, x3 = e_mat.T[ 0 ]
print x1, x2, x3

e_valueは固有値、e_matには固有ベクトルが入っています。e_matに格納されている値がまさにx1, x2, x3となります。それぞれの値は[-0.521333010475 -0.44553760869 -0.727810505167]となっています。

ここで先程の制約条件を思い出してください。平均が0、分散が1ということで定義しました。この条件が合っていることが確認できるかと思います。(厳密には要素数が3ですので、x1, x2, x3にsqrt(3)をそれぞれ掛けた値で分散が1になるのですが、スケーリングするということだけですので気にしなくて構いません。)

もし直線でなく平面でプロットしたい場合は固有ベクトルの1行目と2行目を使います。今一行目でx1, x2, x3としましたが、2行目をy1, y2, y3としてプロットすればよいだけです。3次元の際は、3行目をzとして利用します。

参考文献

[1] http://www.page.sannet.ne.jp/yo-skmt/MultiEx/No13-15MU-Text02.pdf

[2]https://www.nikkei-r.co.jp/glossary/id=1605

おすすめ記事

単純ベイズ分類器(Naive Bayes Classifier )

はじめに

ベイズ定理を利用した分類機を、ナイーブベイズ分類器とよびます。スパムフィルターとして使われており、聞いたことがある人も少なくないと思います。ベイズ分類器は実装が簡単な上に、出力値が確率であるため扱いやすく、ポピュラーな学習モデルです。 当サイトでは長い理論説明はせず、必要最低限の数式をで説明していこうと思います。

具体的なベイズ定理

男性のタバコを吸う人が肺癌になる確率 を求めます。タバコを吸ったら肺がんになるという確率は普通はわかり得ません。ただこの確率が、世の中の肺がん率、喫煙率、肺癌患者の喫煙率が解れば、ベイズ推定で求めることができるのです。

ベイズ式は以下のように表されます。

Pr( A | B ) = Pr( B | A ) * Pr( A ) / Pr( B )

それぞれは次のとおり定義します。

  1. Pr( A ) = 肺癌率
  2. Pr( B ) = 喫煙率
  3. Pr( A | B ) = 喫煙者の人が肺癌になる確率
  4. Pr( B | A ) = 肺癌の人が喫煙者である確率

なお、サイトに寄っては次のようにも書かれたりしています。

P(A|B) =\frac{P(A \cap B )}{P(B)} = \frac{P(B|A)P(A)}{P(B)}

罹患率等のデータの取得

実際に喫煙者の人が肺癌になる確率を求めます。喫煙率、肺がん率、肺癌患者の喫煙率を入手する必要があり、ネットで検索してみました。

男性の肺癌患者の30%はタバコによるガンだそうです[1]。そこで、式4. 喫煙者の人が肺癌になる確率は次のように定義できます。

4. Pr( B | A )=0.3 (喫煙者の人が肺癌になる確率)

続いて、肺癌率は男性の場合は7.4%[2], 男性の喫煙率は29.4%[3]なので, Pr(A), Pr(B)はそれぞれ次のように定義できます。

1. Pr(A) = 0.074(男性肺癌率)

2. Pr( B ) = 0.294 ( 男性喫煙率)

必要な確率は全て知り得たので、最終的なタバコを吸うと肺癌になる確率を求めます。最初の数式にそれぞれの確率を当てはめます。

Pr( A | B ) = 0.3 * 0.074 / 0.294 = 0.075

上記の通り、約7.5%の確率でタバコを吸うと肺癌になることが判明しました。通常、タバコを吸って肺癌になる確率は解るはずはありませんが、ベイズの定理を用いて求めることができるようになるのです。しかしながら、喫煙は8%弱肺癌に寄与するというのはなんとも面白い結果となりました。

ベイズ分類器

ベイズの定理を見たところで、実際にをどうやって分類器として力を発揮できるかを考えます。

ベイズ分類器では、左側のベイズ理論を右側のように仮定していきます。

P(Y|X) = \frac{P(X|Y) \times P(Y)}{P(X)}

上記式は次のようにベイズ分類機では意味すると当てはめています。

Posterior = \frac{Likelifood \times prior }{ evidence }

なお、Likefoodは尤度、Priorは事前分布、Evidenceは周辺尤度と呼ばれています。

尤度であるP(X|Y) 、それから事前分布であるP(Y)は学習データから入手します。今次のような学習データがあったとします。次の学習データセットはエジンバラ大学の教材[4]を参考にしました。以下の表は頻度表といわれます。

天気気温湿度風速ゴルフ
晴れ暑い高いなしNO
晴れ暑い高いありNO
曇り暑い高いなしYES
雨天普通高いなしYES
雨天寒い普通なしYES
雨天寒い普通ありNO
曇り寒い普通ありYES
晴れ普通高いなしNO
晴れ寒い普通なしYES
雨天普通普通なしYES
晴れ普通普通ありYES
曇り普通高いありYES
曇り暑い普通なしYES
雨天普通高いありNO

ここで、X=( 天気、気温、湿度、風速), そしてY=( ゴルフ )となります。YのYESはゴルフした、Noはゴルフをしなかったという意味です

P( X | Y ) と P( Y )を上記の表から当てはめます。

P(Y=YES)=\frac{ YES 数}{ YES 数+NO数} = \frac{9}{14}

P(X=晴れ|Y=Play) = \frac{ YES 数 \cap 晴れ数}{ YES 数} = \frac{2}{9}

P(X|Y)のXは1つだけで説明しましたが、Xが複数条件の時も知りたくなるかもしれません例えば外が晴れ&寒い時にはゴルフするのかどうかなどです。一件わかりますが、寒い&晴れているなどのXの条件数が多くなるほど計算は複雑になります。状態数が多くなるためです。

そこでXの要素はすべて独立しているという仮定を置きます。そうすることで次のように書き換えることが可能となります。

P(X=晴れ,X=寒い| Y ) = P(X=晴れ|Y) P(X=寒い|Y)

独立という仮定をおくことで、分割することができるようになるのです。こうすることでパラメータ数は非常にが少なくなり(具体的には指数関数が線形になり)高次元のデータでも扱うことができるようになります。

データが独立という仮定を置かれている所以がナイーブ(雑)と言われている理由です。

迷惑メールフィルター実装例

実際にベイズ推定を用いた迷惑メールフィルターがどのように作られるのか説明します。

迷惑メール 30通中を受信したとします。迷惑メールの中で次の単語「バイアグラ」、「限定品」、「機械学習」という単語が含まれた表、頻度表を作成します。次のようになりました。

バイアグラ限定品機械学習
回数7170

迷惑メールが30通ですので、ここから尤度表を作成することができます。単に母数で除算するだけです。

バイアグラ限定品機械学習
YES7/3017/300/30
NO23/3013/3030/30

尤度表から、例えば「バイアグラ」が入り、「限定品」が入っていない場合の迷惑メールの確率 P( 迷惑メール | バイアグラ ∧~限定品)を求めることができます。

\frac{P( バイアグラ \cap\lnot限定品 | 迷惑メール) P(迷惑メール)}{ P( バイアグラ\cap \lnot 限定品 )}

さて、ここで∧の条件は独立である場合分解できるということを説明しました。そのため、上記式は次のように分解できます。

\frac{P( バイアグラ| 迷惑メール) P( \lnot 限定品 | 迷惑メール) P(迷惑メール)}{ P( バイアグラ ) P( \lnot 限定品 ) }

続いて迷惑メールでない 正常なメール100通のうち、先程の単語についての尤度表を作成してみます。

バイアグラ限定品機械学習
YES1/1001/1007/100
NO99/10029/10023/100

ゆう土俵から、正常なメールの確率は次のとおりとなります。

\frac{P( バイアグラ| 正常メール) P( \lnot 限定品 | 正常メール) P( 正常 メール)}{ P( バイアグラ ) P( \lnot 限定品 ) }

これで準備が整いました。メールに「バイアグラ」& NOT「限定品」 の条件 で迷惑メールと正常メールのどちらに分類されるかの確率を求めます。今、分母が共通していますので、分子の部分だけに着目して計算をします。(何故分母は無視したのかは後述します)

迷惑メールの尤度= P( バイアグラ| 迷惑メール) P( \lnot 限定品 | 迷惑メール) P(迷惑メール) = 7/30 * 13/30 * 30/130= 0.023

正常メールの尤度= P( バイアグラ| 正常 メール) P( \lnot 限定品 | 正常 メール) P( 正常 メール) = 1/100 * 29/100 * 100/130 = 0.0022

実はこれだけで、この種のメッセージが迷惑メールである確率をもとめることができます。次式で求まります。

0.023/(0.023 + 0.0022 )=0.912

よく見るとわかりますが、算出された合計の尤度で、該当尤度を除算するだけです。結果として、バイアグラを含み、限定品を含まないメールは、おおよそ91.2%が迷惑メールであるということがわかりました。このように尤度をすべて足しこんだ値を、該当の尤度で割ることで、それぞれのメールの確率がわかります。

分母を省略した理由ですが、分母が共通しています。確率を計算するときに打ち消し合って消えるので無視できるのです。

参考文献

  1. https://www.sankei.com/life/news/180104/lif1801040007-n1.html
  2. https://www.haigan.gr.jp/guideline/2018/1/1/180101000000.html
  3. https://ganjoho.jp/reg_stat/statistics/stat/smoking.html
  4. http://www.inf.ed.ac.uk/teaching/courses/inf2b/learnSlides/inf2b12-learnlec06.pdf

次に読むのおすすめ記事

正準相関分析(Canonical Correlation Analysis:CCA)

はじめに

1つの目的変数を、複数の変数から説明するのが重回帰分析でした。正準相関分析は、複数の目的変数に対して、複数の変数がどう説明するのか見る方法です。

重回帰分析の際には、1つの変数を説明するために複数の変数を使っていましたが、正準相関分析は複数の変数を説明するために複数の変数を使います。重回帰分析の一般化、としてみることが出来ます。

数学的な意味の正準相関分析

正準相関分析とは、例えば2つの次元に貼っているベクトル( Va, Ua ), (Vb, Ub )が合った際、同じ空間に射影し、そのベクトルの相関を最大(コサイン積を最大)にするという見方もできる。上図でいえばWa, 3次元から2次元への、Wbは2次元から2次元へのそれぞれ射影行列となる。索引:https://gregorygundersen.com/blog/2018/07/17/cca/

数学正準相関分析は、「複数の変数を2つの合成変数にまとめ、その合成変数間の相関係数を最大にする」という手法です。

つまり数に重み係数をつけて足し合わせたものを変量(正準変量、合成変量とも呼ばれる)として、その変量同士の相関関係である正準相関係数を最大にするような重み相関係数を求めます。

変量は回帰分析においてのy=axにおいてのyの部分です。回帰分析の場合、yの値が決まっていますが、正準相関分析ではyの値は決まっておらず、解析結果によって求めます。

正準相関分析は共通の空間(正準空間)に写像する手法です。

複数の情報源に共通して含まれる情報だけを抽出することによって、情報を統合するための次元圧縮法とも考えることが出来ます。

正準相関分析の具体例

理論だと少しむずかしいので具体的な例で説明します。「数学点数日本史点数握力遠投距離 」 のデータが5人分あったとします。

people1=[90, 40, 50,30]
people2=[80, 50, 55,10]
people3=[30, 80, 40,5 ]
people4=[50, 50, 35,50]
people5=[70, 80, 30,40]

同じ人たちの「走り幅跳び距離、英語点数」のデータもあります。

people1=[1.5, 40]
people2=[3.1, 30]
people3=[1.7, 50]
people4=[0.5, 80]
people5=[1.3, 90]

走り幅跳び距離と英語点数が、どう影響を受けるかということを考えることができます。

2つグループの相関関係を無理矢理最大化するような式をそれぞれ作ってやる、という事が正準相関分析となります。先に説明したとおり、解析結果後にそれぞれのグループで合成変量がでてきます。y=aXでいうところのyの部分です。このyの値をどう捉えるかはその観測者次第となります。

分析結果については、PCAのように寄与率と、最も相関があるように出来た重みが出てきます。N次元のデータであれば、N個の軸での解析が可能です。ただ、2つのグループの次元が違う場合がほとんどです。その場合は、最も長い次元にあわされますが、経験的に5−6個の軸から見れば大体解析できます。以下にCCAのコードを記載します。(察しの通り、正方行列ではないので特異値分解を中でおこなっています)

from sklearn.cross_decomposition import CCA
    cca = CCA(n_components=5)
    cca.fit(X1, X2) 

n_componentsは、寄与率の高い上から何個までの重み群と寄与率取ってくるのか、指定するものです。解析した異次元数を表しています。

CCAが中で何をしているかは、数学的な説明が多くなる為詳細を省きます。正準相関分析に関して解き方などの理論が気になる方は、君山 由良先生の多変量回帰分析・正準相関分析・多変量分散分析―多変量間の相関と因果関係の因子 (統計解説書シリーズ A-16) をおすすめします。非常に丁寧に解説されております。

参考文献

おすすめサイト

回帰分析(重回帰分析)

回帰分析の概念はとても簡単です。データにフィッティングする直線を見つけるという事だけです。例えば、身長と体重の相関関係を知りたいとします。直感的にですが、身長が高くになるに連れて、体重も増えそうですね。

体重 = a * 身長

となるような直線を推定する、これが回帰分析です。ここでよく見ると、y = a * xという式と同じです。( yが体重、xが身長)。つまり、回帰分析は傾きaを求めるための分析方法ということです。

さて、今の回帰分析は単回帰と呼ばれています。もっと複雑に体重を求めたいかもしれません。例えば説明変数「性別、宗教、胸囲、頭の多さ、身長」から、「体重」を求める、というものです。このように複数の説明変数から、体重を求めたい場合(説明変数の傾きを求めたい場合)は重回帰分析と呼ばれています。式は次のような形となります。

y = a_1x_1 + a_2x_2 + a_3x_3 + ...

説明の通り複数の説明変数がある場合が、重回帰分析と呼ばれています。重回帰分析でどうやって係数群をもとめるかは最小二乗法のコラムをご覧ください。

少し話はそれますが先程上げた「性別、宗教、胸囲、頭の多さ、身長」から「体重」を求めるという行為は実はやってはいけません。コンプライアンス的にという意味ではなく、重回帰分析のルールとして次のがあるためです。重回帰分析には次のルールがあります。

  • 説明変数は独立であること。
  • 説明変数は連続値であること。

身長と、胸囲は相関関係があり独立ではありませんね。身長が高くなると胸囲があがるという相関関係がありそうです。また、性別は0か1ですので連続値ではありません。宗教もイスラム、キリスト、仏教、神など連続ではないです。

なお、ルールを無視しても精度は出ます。ただ適切ではないです。論文などの公式なところに使えば突っ込まれるだけで、社内のデータなどでは使っても問題はありません。それでは、多変量解析のデータを解析する上でどういう方法を使えばいいかについて、迷った場合は以下を参照ください。

従属変数従属変数が量的独立変数が量的使用できる多変量解析の種類
ありYESYES重回帰分析
ありYESNO(質的)分散分析、多重分類分析、数量化理論Ⅰ類
ありNOYES重回帰分析、判別分析、ロジスティック回帰分析
ありNONO(質的)対数線形モデル、数量化理論Ⅱ類
なしYES因子分析、主成分分析
なしNO(質的)クラスター分析、多次元尺度法、数量化理論Ⅲ類

従属変数とは、y=a_1x_1 a_2x_2でいうところの、yが存在するかどうかを表します。質的なデータとは0,1で区別されるようなデータ、量的とは連続なデータか、0〜1の範囲でくまなく分布する連続するデータか、という意味を示します。

例えば、宗教(質的)と性別(質的)から癌になる確率(量的)を求めたい場合は数量化理論Ⅰ類を用いることとなります。

さて、重回帰分析では、一個の従属変数と複数個の独立変数の線形合成変数の相関が最大となるような重みを求めることを目的としています。(従属変数y、独立変数x、重みa)説明変数が一つなら単回帰分析。2つ以上なら重回帰分析となり、単回帰分析とは直線近似を指し、重回帰分析はそれの一般化で超平面を指します。

今、目的変数yに対し、説明変数xがどれだけの重みを持っているのかということを求めるのが重回帰、単回帰分析であると説明しましたが、改めて式にして表せば、

(単回帰分析) Y=aX

(重回帰分析) Y=aX // 太字になっているということで、ベクトルを表している。

終わりに

目的変数がひとつであることが回帰分析であると説明しました。重回帰分析は単回帰分析の一般化です。重回帰分析を更に一般化すると、正準相関分析となります。正準相関分析については、また別途解説いたします。

参考文献

https://ja.wikipedia.org/wiki/%E5%9B%9E%E5%B8%B0%E5%88%86%E6%9E%90

おすすめの記事

ロジスティック回帰分析

ロジスティック回帰分析とは、モデルがlog( \frac{1} {1 - p} ) = a_1 x_1 + a_2 x_2 + a_3 x_3 ...となるような{a_1,a_2,a_3}ベクトルを探すモデルです。(上記pはベクトル、x1..xnは今既知の値です。)

一旦モデルを作成さえすれば、あとは実測値を入れるだけで良いのがポイントで、このあたりは回帰分析と同じとなります。それでは一体何が違うのでしょうか。

ロジスティック回帰分析では、2値のデータの予測をしたいときに利用します。(男か女か、イスラム教かそれ以外か、日本人か外人かなど)。回帰分析では、値の推定に正規分布を利用しているため、正規分布に従わない、たとえば0か1かしか取らないものの推定が出来ないのです。少し専門的な言い方をすると、量的変数の解析には回帰分析、質的変数の解析(ただし2値に限る)にはロジスティック回帰分析を使います。量的変数は連続値のデータで、質的変数は離散的なデータと捉えて問題ありません。

出力が2値ですので、面白い使い方として、たとえばこれをしたら売れるか(1)、売れないか(0)みたいな出力結果を得られます。ですので、例えばマーケッターの方があるターゲットに売れるか売れないかの2値をとる計算などにロジスティック回帰分析が使えます。

以下のような関数を使います。

P = exp( Z ) / ( 1 + exp( Z ) ) =  1 / ( 1 + exp( Z ) )

ここで、上記の中のZは次のように表します。

Z = B_0 + B_1*x_1+ B_2*x_2+ B_3*x_3 …

補足ですがPの式はよくみるとシグモイド関数です。このBの係数(偏回帰係数
とよばれる)を求めるために式を変形します。

ln( P / (1-P )) = Z

この式をみるとわかるとおり、あとは、最小二乗法を用いて偏回帰係数を求めるだけで、回帰分析と全く一緒となります。これにより、ロジスティック回帰分析が出来ます。なお、最尤法(コンピューターをぶん回してフィッティング関数を見つける方法)を利用する方法でも解は求まるのですが、ほとんど最小2乗法を利用した計算法で事が足ります。

参考文献

おすすめの記事

最尤法

最尤法は、ある観察された離散的なデータを、連続的な関数にフィッティング(推定)する事を目的とした方法です。

今観測されたデータを正規分布や二項分布、ガンマ、ポアソンなど何でもいいのですが、連続関数にうまく当てはめてやる。うまく当てはまったということは、その際にその分布のパラメータの推定を行えていることを意味します。

今自分で適当に仮定した確率モデル(例えば正規分布など)があった場合、観測された実際の分布との適合度を示すものを「尤度」と呼びます。

尤度は高ければ高いほど、その仮定した確率モデルのフィッティングがうまく言っている事を表す。(尤度は、ある理論分布(正規分布など)を仮定したとき、現実に得られた実データと 同じ結果がその分布から生じる確率)

仮定する確率モデルのことを「尤度関数」と呼ぶ。(つまり観測結果をフィッティングした関数のこと)

尤度関数は単純に正規分布の場合や、複数の分布の混合の形で表される場合もあります。

最尤推定、最尤法とは、尤度を最大化するような尤度関数を決めてやる(パラメータを推定する)ことです。

参考文献

https://ja.wikipedia.org/wiki/%E6%9C%80%E5%B0%A4%E6%8E%A8%E5%AE%9A

おすすめの記事

FPマイニング(頻出パターンマイニング )

FPマイニングとはFrequency Pattern Mining の略で、物事の関連を見出すことを目的とします。たとえば、「おむつを買った人とはビールを買う」などそうした関連を見出すことを目的とします。どういうところで活躍しているのでしょうか?たとえば広告の技術などで活躍します。例えば、AというサイトをみてBというサイトを見たい人はCという商品が欲しい。それではCの広告を出そう!という事になります。こうした関連性を見出したいときはFPマイニングが利用されます。それでは説明に入っていきます。

アソシエーションルール

「アソシエーションルール」という言葉は、日本語で、連関ルール、関連ルール、相関ルール(⇐ 全部同じ意味を表す)と呼ばれています。記号で表すと以下のようになります。

{X}→{Y}

XであればYであるという意味で。よく例題で用いられるのはビールと紙おむつです。

  {ビール} →{紙おむつ}

上記の意味は「ビールを買った人はおむつを買う」という意味です。

FPマイニングはそのまんま実装すれば非常に簡単ですが問題があります。データ量が多くなると、とたんに難しくなるということです。これは、全ての組み合わせを探索する必要があるため、組み合わせ爆発が発生することによります。(要素数をXとすると2Xパターンが存在)

そのため工夫としてAprioriアルゴリズムやFP-Growth といった手法がとられており、実際にはFP-Treeを構築するFPGrowthがというアルゴリズムが用いられています。

アソシエーションルールの評価方法について

アソシエーションルールの評価方法には全部で3つの指標があります。

  1. 支持度( support )
  2. 確信度( confidence )
  3. リフト( lift )

結局は確信度が大事なわけなのですが、確信度は支持度(Support, sup)で算出されるため、Supを効率的に算出することを目標としています。

まずは、上記1,2,3それぞれの式がどのように表されるか説明いたします。

1. ルールX ⇒ Y の支持度

supp(X ⇒ Y) = σ(X ∪ Y)/ M

2. ルールX ⇒ Y の確信度

conf (X ⇒ Y) = σ(X ∪ Y)/ σ(X) もしくは、

conf (X ⇒ Y) = supp(X ⇒ Y) / supp(X)

3. ルールX ⇒ Y のリフト

lift(X ⇒ Y) = conf (X ⇒ Y) / supp(Y)

ここで、上記のσの意味は全体の集合から、要素Xを含むのが幾つあったのかを表すために用いられています。

例えば、集合が3つあったとして、それぞれの要素X, Y, Z, Wは次のようになっているとします。

{X, Y}, {X,Z}, {X, Y, Z, W} 

ここで、σ( X )は 3となります。そして、XとYをふくんでいる集合の数はσ( X ∪ Y ) = 2 となります。なお、集合の数は3つなので、式中にあるMは3となります。
それぞれの「支持」「確信」「リフト」に関して説明は以下のとおりです。

  • 支持度については全体の集合からそのルールが幾つあったのかという確率である。
  • 確信度については、ルールXがあったときに Yも発生する確率である。(XであればYも買う確率)
  • リフトについては 確信度を支持度で割っている。支持度がひくければ低いほど値が高く(データの中から出現度合いが低いとリフトが上がる)、確信も高いと高くなる。これは、確率P r ( X∪Y ) /Pr(X)Pr(Y)の近似値としても知られている。ちなみ、これは相関である。つまりLIFTは相関を表している
    • Liftは同じ確信があるときに、どちらの商品が人気あがあるのかがよくわかる。例えば、0.3/0.1, 0.3 / 0.05だと後者の方が大きい。除数は支持率(全体の数である)

3つの指標において、最終的に関連が高いものは次の二つの事が言えます。

  • 確信度conf が高い
  • 支持度support が高い

FPマイニングはこうした数値を出して、確信度が高い!支持度が低かった!など分析していくことです。

前述のとおり、全探索すれば関連性は見いだせるのですが全てでるので良いのですが、実際には不可能です。そのため、以下のような簡略化するアルゴリズムが用いられます。

Support(頻度)さえでれば、Confidenceが求まることになります。Confidence( ビール →おむつ ) = support( {ビール、おむつ} ) / support( ビール )

ですので、どのように効率的にSupportをだすかが鍵となります。

Apriori アルゴリズム

全データ数を見て、その関連を見つけるのは大変です。(組み合わせ爆発問題)。そのためAprioriアルゴリズムによって効率的に上記のアソシエーションルールを満たすものを探し出します。以下はその様子を図示した内容です。

集合から、要素を取り出し、その数をカウントしてC1を作成します。しきい値(この場合は1とした)より低いものを無視して、L1を作ります。L1の要素のすべての組み合わせからC2を作ります。ここで、C2のセットで1つ減らした要素( {1 2}であれば{1}と{2}が) L1に入っていない場合、C2の要素としてカウントは出来ません。C2作成に当たりその例はありませんでした。ただC3でその事例が発生します。

C2からスキャンして、先程と同様にしきい値より下のものを削除しL2をつくります。以下同様にC3をつくります。先に説明したとおり、C3の部分集合がL2に含まれていない場合は、その部分集合は入れてはいけません。そのため、{2 3 5}のみが候補として残ります。

さて、Supportさえ出てしまえば、Confidenceを算出することが出来ます。Aprioriアルゴリズムの問題点はメモリです。100個のアイテムから候補を作成するのに2の100乗の候補が必要になり、メモリを圧迫するうえ、データのスキャンにも時間が必要となります。

参考文献

おすすめの記事

主成分分析(PCA)

コンピューターサイエンスで非常に重要な技術の一つです。次元圧縮の定番で「まずはPCAでもしてみるか」的な感じで利用されます。PCAとはPricipal Component Analysisと呼ばれ日本語では主成分分析と呼ばれます。

コンピュータサイエンスを専攻する学生にとっては、必ず学んでおいた方がいい技術の1つです。ただ殆どのサイトは数式を並べていて、結局何なのかがわかりづらいな、と思います。機械学習、統計学にしろ概念が非常に重要です。このブログでは概念を主に説明します。最後に簡単な実装例をおつけしておきます。

PCAとは(概念)

PCAは次元圧縮や特徴抽出というふうに言われます。なぜでしょうか。

それは、たくさんのデータから、そのデータを構成する基礎となるようなデータ(軸)が求まるからです。

少し詳細に説明します。例えば、ものすごい次元のデータが膨大にあったとします。
全部保持するのは容量的に無理です。ここで基礎となるような数個のデータと係数で、その莫大な集合をほぼ表現できるとしたら便利だし、これは次元圧縮と言えるのではないでしょうか。そして基礎となるような数個のデータが膨大なデータの基礎となるなら、これはデータの特徴とも呼べそうです。PCAはそれをするため次元圧縮や特徴抽出と言われています。

基礎となるようなデータは主成分ベクトルとよばれます。主成分ベクトルが面白いのは基底ベクトルであることです。基底ベクトルは、互いに独立で直交しているという特性をもっているベクトルで、その座標系での軸となるベクトルです。

実はPCAにより算出される主成分ベクトルはすべて基底ベクトルです。PCAにより算出された基底ベクトルが織り成す特徴量空間は元データがよく見える空間となります。(分散してみえるようになる)。こうした基となる軸のベクトル(データ)がPCAにより求まるのです。

ちなみに、基底ベクトルであれば掛け合わせるだけで簡単にその基底ベクトルを軸とした空間に変換が出来ます。線形代数の話になりますが、直交する基底ベクトルのみで構成されている行列を直交行列といい、その直交行列をデータにかけると、直交行列を系とした座標系に変換ができます。

PCAにより得られる主成分ベクトル同士は直交しており、それぞれの軸からでも分散をが大きくみえますが、その量は軸別で分散値が異なります。このどれくらい分散しているか(よく見えるか)というのは、寄与率で判断ができます。寄与率が高いほど主成分ベクトルがどれくらい重要かが判定できます。

蛇足で、別に直交していなくてもよく見える軸で見ればいいではないのか、というようなアイディアがあり、これはICA(独立成分分析)とよばれ音声処理などで活躍しています。

PCAの具体的な処理手順の

主成分分析とは、以下の手順をするだけで求まります。

・共分散行列を求めて固有値問題を解く

これがPCAです。従って、問題は以下の2つになります。

  1. 共分散行列を求める。
  2. 固有値問題を解く。

この二つのプロセスで解けるということで、わかってしまえば非常に簡単です。

共分散行列とは

共分散行列は、分散共分散行列とも呼ばれています。共分散行列の求め方については以下の具体例で示します。今、N次元のV1, V2, V3, V4のベクトルがあるとするとします。

そしてV1,V2,V3,V4の平均ベクトルV_avrを求めます。

V_{avr} =\frac{ V_1 + V_2 + V_3 + V_4 }{ 4}

続いて、V1, V2, V3, V4それぞれの共分散行列を作ります。例えば、V1の分散行列は次のようにして作成されます。

(V_1 - V_{avr}) \times (V_1 - V_{avr})^T ( ^Tは転置行列を表します。結果的にNxNのマトリックスが作成されます。)

例えば、V1が[ 2, 3, 4]という行列で平均ベクトルが[1, 1, 1 ] だったとすれば、V_1 - V_{avr}を行い[1,2,3]というベクトルを得ます。後は式にならって、


{MAT}_1 = \begin{bmatrix}1 & 2 & 3 \\ 2& 4&6\\3&6&9\end{bmatrix}

となり、V1の共分散行列MAT1が求まりました。V2, V3, V4についても同様のことを行います。

最終的な共分散行列をもとめるには、それぞれの共分散行列を全て足し、ベクトルの数(4つ)で割ればそれで終わりです。

共分散行列=(MAT1 + MAT2 + MAT3 + MAT4)/4 

これで目的の共分散行列を作ることができました。意外と簡単です^^;

固有値問題とは

固有値問題は以下のようなベクトルと値を求める操作です。

Ax = λx

ここで、λは固有値とxは固有ベクトルと呼ばれます。ここでAは行列を表します。 ここでいうAが先程のNxNの行列(共分散行列)となります。

固有値問題を解く方法はこの範囲を超えますので、様々なサイトをご覧ください。

PythonのLinear Algebraのライブラリを使って解く擬似コードを紹介します。

import numpy
from numpy import array
from numpy import linalg
def PCA(X):
  # 共分散行列の計算
    x_mean = array( map( sum, X.T) ) / len( X )
    X_ = X - x_mean  
    X_t = numpy.dot( X_.T, X_ ) / len(X)
    #固有値、固有分解
    lam, vec = linalg.eig( X_t ) # <== 固有値、固有分解が求まる!!
    ans = zip( lam, vec.T )
    ans.sort( reverse = True ) # <= 固有値が高い順にソート

linalg.eigにより固有値分解、固有値ベクトルを算出します。固有値分解の方法については興味があれば自分で調べてみてください。このように、固有値分解のライブラリさえ手に入っている状態であれば簡単にPCAを行うことができます。

ランク落ち

固有値分解を行う上で一つ注意がRANK落ちに注意してください。ランクがFULLでないと(つまり4x4の行列であれば、4つの独立したベクトルにより、その行列が構成されていないと)不正な値が出ます。(使うライブラリや実装方法によりますが)

そうした場合、正則化させてランク落ちをなくすことや、行列のランクがFULLでないのであれば特異値分解により求める必要があります。

参考文献

https://ja.wikipedia.org/wiki/%E4%B8%BB%E6%88%90%E5%88%86%E5%88%86%E6%9E%90

おすすめ記事