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

おすすめ記事

ページランク( Page Rank )

はじめに

ページランクの仕組み 、アルゴリズムを意外と知らない人は多いと思います。ページランクアルゴリズムはGoogleの創業者であるLawrence Page とSergey Brinにより考案されました。今や世界を支えていると言われる検索エンジンアルゴリズムはとてもシンプルな式でした。

PR(webPage_x) = (1-d) + d \times (\frac{ PR(webPage_1)}{Count(webPage_1)} + ... + \frac{ PR(webPage_n)}{ Count(webPage_n)})

PR()はページランクを表します。最初は定まっていないので初期値として全てのページに適当な値を入れますが1を入れるケースが多いです。

dは減衰率でオリジナルの論文では0.85として設定されています。

Count()はWebページが持つリンクの数です。たとえばWebPage1がリンクを10個持っていればC( WebPage1 ) は10となります。

上記式をみると、自分のページランクをリンク数で除算し分け与えているという意味となります。こんな単純なアルゴリズムから世界を変えるほどの技術となっているのは本当に驚きです。

さて実際の計算しページランクがどのように更新されるのか見てみます。

Page Rankを計算してみる

以下のようなページがあったと仮定します。

AはB,C,Dへバックリンクを貼る、BはDへ、CはA, Dへ。DはA,Cにリンクを貼る。それぞれ、Count(A) = 3, Count( B ) = 1, Count( C ) = 2, Count( D ) = 2となる。

今、4つのウェブページがあり、それぞれリンクを張っています。Aは3つ、Bは1つ、Cは2つ、Dは2つリンクを張っています。このリンクを貼られたという関係を表すと次のようになります

AからB から C から D から
Aへ0011
Bへ 1000
Cへ 1001
Dへ 1110

例えば、Aは(B,C,D)にリンクを張っていますので(B,C,D)にそれぞれ1を立てました。

続いて確率モデルの表を先の表から作成します。今Aのリンク数、Count(A)は3です。Aが他のサイト(B,C,D)への分配数はその逆数である1/3ずつが分配されます。上の表をそれぞれのリンク数で除算し,作成された表は次のとおりです。

AからB から C から D から
Aへ00\frac{1}{2}\frac{1}{2}
Bへ \frac{1}{3}000
Cへ \frac{1}{3}00\frac{1}{2}
Dへ \frac{1}{3}1\frac{1}{2}0

行列の縦軸に注目すると、合計が必ず1になっていることに気づくと思います。また、マイナスの値が入っておりません。一般的にこうした行列のことを確率行列( column-stochastic matrix)と呼びます。

ページランクベクトルを用意し、先の式に適用していくとページランクが出てきます。今、初期値のページランクを1とします。

PR(A) = 1, PR(B) = 1, PR(C)=1, PR(D)=1

準備がととのったのでページランクを計算します。わかりやすさを重視するために愚直なコード、いわゆるベタ書きコードで表します。実際のプログラムでは、100万ページがあってもエレガントに解けるようなソースコードの書き方を目指してください。

def main():
    # PageRank初期値                                                                                                                                                                                                              
    PR = {'A': 1.0, 'B': 1.0, 'C':1.0, 'D': 1.0 }
    # リンク数の逆数(確率)                                                                                                                                                                                                      
    Prob = {'A': 1/3.0, 'B': 1, 'C':1/2.0, 'D': 1/2.0 }
    # 減衰率                                                                                                                                                                                                                      
    d = 0.85
    # 10回回す                                                                                                                                                                                                                    
    for i in range( 10 ):
        new_PR_A = PR[ 'C' ] * Prob[ 'C' ] + PR[ 'D' ] * Prob[ 'D' ]
        new_PR_B = PR[ 'A' ] * Prob[ 'A' ]
        new_PR_C = PR[ 'A' ] * Prob[ 'A' ] + PR[ 'D' ] * Prob[ 'D' ]
        new_PR_D = PR[ 'A' ] * Prob[ 'A' ] + PR[ 'B' ] * Prob[ 'B' ] + PR[ 'C' ] * Prob[ 'C' ]

        PR[ 'A' ] = (1-d) + d * new_PR_A
        PR[ 'B' ] = (1-d) + d * new_PR_B
        PR[ 'C' ] = (1-d) + d * new_PR_C
        PR[ 'D' ] = (1-d) + d * new_PR_D
    # 新しい値を出力
        print( "%.4f %.4f %.4f %.4f" % (PR[ 'A' ], PR[ 'B' ], PR[ 'C' ], PR[ 'D' ] ) )
if __name__ == '__main__':
    main()

10回ぶん回した結果は次のとおりです。段々と値が安定しているのがわかります。

ページ数や構造がシンプルな例でページランクの計算方法をみてみましたが、実装がシンプルであり、とてもわかり易いアルゴリズムであると思います。ただ、世の中のモデルがこのように単純な問題であるとは限りません。例えばリンクを貼られていないページがあるとどのような問題がでてくるのでしょうか。

Rank Sinks問題

あるページがリンクを全く持っていない問題のことをRank sinks 問題と呼びます。例えば次のようなケースの場合です。

Dが最後になっている場合、Rank Sinksと言われています。流しのシンクのようになっているからです。Dはリンクを持っていないので移動する確率は0であり、周りからページランクをもらうだけです。Dは必然とページランクが上がる一方で、他に配る宛もなく、これは公平ではありません。

この場合はシンクDの移動確率は、全ページ数の逆数を利用することで問題解決を図ります。今(A,B,C,D)と総ページ数4ですので1/4がDの移動確率として定義されます。

AからB から C から D から
Aへ00\frac{1}{2}\frac{1}{4}
Bへ \frac{1}{3}00\frac{1}{4}
Cへ \frac{1}{3}00\frac{1}{4}
Dへ \frac{1}{3}1\frac{1}{2} \frac{1}{4}

全ページ数の逆数である理由は、最後にたどり着いたページからランダムに別のページに行く確率は等確率である、という仮定に基づいています。

注意(ページランクのもう一つの式

ページランクの式としてはNで割る場合もあります。

PR(A) = \frac{(1-d)}{N} + d \times (\frac{ PR(webPage_1)}{Count(webPage_1)} + ... + \frac{ PR(webPage_n)}{ Count(webPage_n)})

最初の項目でNという項目が入っただけです。Nはページ数を表しており、先の例であれば4が入ります。この修正された式の理由は全てのページランクの合計は1である、つまり確率モデルとして正しい形を求めた為です。従って本質的な意味の違いはありません

ところで、PageRankのアルゴリズムはマルコフ連鎖であるといわます。マルコフ連鎖とはざっくりいうと、次の状態は前の状態からのみ決定されるというものです。ページランクのアルゴリズムを観ると分かる通り、一つ前の状態からしか次のページランクは影響をうけません。そのため、マルコフ連鎖を仮定してという言い方がされます。

ページランクの仕組みは理解できましたでしょうか?当サイトが皆さんの知識の工場に役立ってもらえれば幸いです。

参考サイト

関連記事

単純ベイズ分類器(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

次に読むのおすすめ記事

Matrix Factorization(MF)

はじめに

Matrix Factorization(MF)について解説していきます。

現在世にある推薦システムは殆どMFをベースに作られているといっても過言ではありません。Simon Funkという人がNetflix Prizeに応募するために開発し、とても話題になりました。SVDを用いて推薦システムを構築したのですが、あるトリックを用いたSVDであったため、中のエンジンはFunk-SVDなんて言われたりしています。他にもNMFを利用したりSVD++など改良が進んでいますが、メインのやり方は本稿の説明でわかるかと思います。

結局、やることは従来のSVDとおなじになる。ただし、求め方で工夫を行う

何故SVDを改良する必要があったのかと言えば、SVDは疎でない(スッカスカでない)行列の解析には素晴らしく良い結果をもたらしていました。しかしながら、商品や映画の推薦となると、全てのユーザーが全てのアイテムを推薦している訳はありません。実際に映画のレーティングの行列を作ると95%のがNaN(未評価)であったということからも、推薦のもとになる行列はとても疎な行列です。この状況を改善しようとした試みがFunkSVDだったというわけです。

概要

Matrix Factorizationの流れを説明します。次のようなレーティングマトリックスを2つ(Item & User )に分けるものです。他のサイトと差別化するために、本ブログでは 一番わかり易い 表ベースの解説をします。今、以下のようなレーティングマトリックスがあったとします。

Movie 1Movie 2Movie 3Movie 4
User 15NaN83
User 277NaN4
User 3NaN677
User 410NaNNaNNaN

上記の行列を2つの行列q, pに分ける事をMFでは行います。

\hat{r}_{ui} = q_i^Tp_u

\hat{r}は推定されたレーティングという意味し、実際のレーティングはrとします。

p_uはユーザー、q_iはアイテムを表す行列です。今回の目的は、ユーザー特徴量行列とアイテム特徴量行列をつくり、それらを掛けると元の行列にする事です。

これは実際のレーティングrと推定レーティング\hat{r}の2乗誤差を最小にするというな行列を分解する,という意味でもあります。式で表せば次のようになります。

minimum( p,q ) \sum_{u,i\in K}(r_{ui}-q_i^Tp_u)^2

p, qを作る

ユーザ特徴量行列p、アイテム特徴量行列qは今は未知なので適当なランダム値で埋めてやります。

p=
Factor 1Factor 2Factor 3
User 11.40.4-1.2
User 21.31.5-0.3
User 30.80.10.6
User 4-0.70.51.1
q=
Movie 1Movie 2Movie 3Movie 4
Factor 10.60.51.3-0.5
Factor 2-1.1-0.40.71.3
Factor 310.30.5-0.1

さて、説明のために、User1のMovie3に注目します。

実際の評価rはレーティングマトリックスをみてみると”8″でした。

ユーザ特徴量行列p(青)ではUser 1は[1.4, 0.4, -1.2]

アイテム特徴量行列q(緑)のMovie3は[1.3, 0.7, 0.5]

です。推定された評価\hat{r}は[1.4 * 1.3 + 0.4 * 0.7 + 0.5 * -1.2 ] = 1.6となります。なお、算出しなくてもいいのですが2乗誤差はは(8 - 1.6)^2=40.96となりました。

更新処理

誤差を元にp, qの値を更新します。p,qの行列から参照した要素は次のとおりでした。

From p1.40.4-1.2
From q1.30.70.5

新しい特徴値は次の計算式を用いて求めます。

p_i = p_i + 2 \alpha ( r - \hat{r} ) q_i

q_i = q_i + 2 \alpha ( r - \hat{r} ) p_i

ここでαは学習率と呼ばれ、0.1くらいにするケースが多いです。これに習ってpの最初の要素を計算してみると、次のようになります。(p_0は1.4, q_0は1.3を指します)。何故この更新式になるかは[2]を参照ください。

1.4 + 2 \times 0.1 \times ( 8 - 1.6 ) \times 1.3 = 3.064

出た値を用いてp_0の値を更新します。

p3.0640.4-1.2

同様にp, qの残りに対しても更新処理を行います。なお、qの更新式ではpの値も使いますが、既に更新された値を使います。

上記例ではp_0の値は3.064になりましたので、q_0を求める際には、新しい数値を使います。あとは繰り返し処理を繰り返し行い、 すべてのp, q行列を更新していきます。以上となります。

これを学習データを用いて適当に取ってきた値で特徴量行列を更新していきます。

ちなみに本当の評価rがNaN値であれば重みの更新処理をしなくても構いません。(というかしません。)。これがNaN値を考慮したSVDであると言える理由です。

重み更新は、ニューラルネットワークなどではおなじみなのですがStochastic Gradient Descent (SGD):確率的勾配降下法と呼ばれています。ランダムに(確率的に)適当な値を引っ張ってきて、重みを更新している(勾配をもとに、2乗誤差を低くするように降下させている)からです。差分がゼロ(微分がゼロ)になれば、最適解に到達したと言えます。

なぜ特異値分解:SVDと呼ばれるのか

今までの操作のどこが特異値分解なのでしょうか。

特異値分解とは、ある行列をエラーを最小とするように2つに分割する操作でした。今回分割する際に行ったことは2乗誤差が最小になるように分割していく操作です。SVDはもともと最小二乗近似を得る方法とも言えます。 実際に開発者のFunkも自身のブログで次のように述べています。

Singular value decomposition is just a mathematical trick for finding those two smaller matrices which minimize the resulting approximation error–specifically the mean squared error (rather convenient!).

さて、大まかな説明はこれで終わりなのですが、更にいろいろな工夫をすることで精度が高くなります。SVDのかわりにNMFを使うことだったり、今回はマトリックスを分割する際のエラーの計算方法として単純にpqに分けるようにしましたが、実際に精度を出すためには正則化したり、バイアス項を突っ込む必要があります。

さらなる改良

正則化

結局重みは次の式を最小化するということでした。

minimum( p,q ) \sum_{u,i\in K}(r_{ui}-q_i^Tp_u)^2

ここで上記式に正則化をします。(過学習を防ぎ、一般的に使えるようにするためにペナルティ項をおく)

minimum( p,q ) \sum_{u,i\in K}(r_{ui}-q_i^Tp_u)^2 + \lambda (\|q_i\|^2 +\|p_u\|^2)

最初の項目に\lambda (\|q_i\|^2 +\|p_u\|^2)をおいただけです。この意味ですが、1つの映画を酷評したユーザーがいて、その映画以外は評価していなかったとします。すると、このユーザーから他の映画へのすべての評価が非常に低くなります。この正則化をおけば、 例えば q_iに大きな値を与えることを抑制することになります。

バイアス項

推測されるレーティングは次の式でもとまるものでした。

\hat{r} = q_i^Tp_u

ここに3つのバイアスを加えます。

  • \mu 全アイテムの平均レーティング
  • b_i アイテムiの平均レーティングから\mu を引いたもの
  • b_u はユーザーuの平均レーティングから\mu を引いたもの

最終的に、レーティングの式は

\hat{r} = q_i^Tp_u+\mu + b_i + b_u

となります。そして、こちらにも正則化項をいれます。

minimum( p,q ) \sum_{u,i\in K}(r_{ui}-q_i^Tp_u)^2 + \lambda (\|q_i\|^2 +\|p_u\|^2+b_i^2+b_u^2)

これが最終的なSVDの式となります。長くなりましたが、これでMFの説明は終わりとなります。

最後に

その他、NMFやSVD++など、エンジンを変えて更新式を実行したい場合は、次のURLに詳細に書かれています。[1]

もしNMFによる実装方法などを記事化してほしい方がいれば、コメントをいただけますと幸いです。

索引

関連記事

特異値分解(SVD)

特異値分解とは

特異値分解はPCAとほとんど同じと思ってください。SVDはPCAが適用できないような行列に対して、データ削減や特徴抽出が可能な方法となります。

PCAを行う行列は正方行列であり、ランクがフルである、つまり逆行列を持つ正則行列である必要がありました。SVDでは疑似逆行列を利用するので正方行列ではない行列やランク落ちしている行列からでも、計算処理が可能となります。

MAT_A = \begin{pmatrix} a & b \\ c & d \\ \end{pmatrix}

正則行列は、逆行列を持つ正方行列を指します。正方行列とは上記のようにN x N行列のことで、横と縦の数が一致した行列を指します。

データは常に正則行列であるとは限りません。PCAは正則行列である必要がある制約条件がありました。SVDは疑似逆行列を利用することで、PCAに似た計算処理を可能にするのです。

どのようにしてSVDを用いて推薦システムを作るのか解説していきます。

SVDによる推薦システム

今ユーザー3名が5つの映画にそれぞれ次のような評価をつけたとします。

review_user1 = [1,4,5,4,5]
review_user2 = [1,3,1,4,4]
review_user3 = [1,2,5,4,1]
mat = numpy.array( [review_user1, review_user2, review_user3 ] )

上記のデータをSVDします。numpyのライブラリを使うことで一発で求まることができます。

def learn( mat ):
    u,s,v=numpy.linalg.svd( mat )
    return [u, s, v]

U, S, Vが求まりました。この関係ですが、元の行列matはU S V^tで求めるように分解をしています。Sは 対角行列 となっており、対角成分に特異値が並んでいます。対角成分は左上が最大で右下が最小となるようになっています。なんだか面白い行列が出てきました。

さて、このUSVを求めたら類似度はどのように求まるかというと、つぎのようにすれば求まります。eには現在自分の映画5個における評価を入れておきます。e= numpy.array([1,5,4,4,5])

def getSimilar( usv, e, dim =2):
    u, s, v = usv
    vv = numpy.array( numpy.matrix(v).T )    
    U = numpy.matrix( u[ :, 0:dim ] )
    V = numpy.matrix( vv[ :, 0:dim ] )
    S = numpy.matrix( numpy.diag( s[:dim] ) )
    E = numpy.matrix( e )
    pos = E * U * S.I
    result = []
    for i, user in enumerate(V):
        val = float(user *pos.T)
        val /= (numpy.linalg.norm( user ) * numpy.linalg.norm( pos ))
        result.append( [val+1, i ])
    return sorted( result, reverse=True )

getSimilar関数はコサイン積をとっています。Dimは次元数という意味で、PCAで言うところの寄与数の高いところ順に内積を取るイメージを取ってください。これにより自分と似ているユーザーが返ってきます。

最後に

実は今紹介したSVDによる推薦システムは、データにSVDを適用し、コサイン積を求めたものでした。推薦システムとしては1つ大きな問題が残っております。それは評価していないアイテム、例えば見ていない映画の評価を0とするとした評価がSVD適用時に考慮されてしまっているということです。見ていないものが、つまらなかった(評価0)は少し辺です。そして、スッカスカの巨大行列をSVDするっていうのは何か違和感があります。そこで評価していないものはスキップするというようなアイディアを取り入れた方法を使う必要があります。これが matrix factorization(MF)という手法になります。Netflixの推薦エンジンのコンペ
で圧倒的なスコアを叩き出したエンジンがまさにMFを使っていました。MFについて、次回ご説明します。

参考文献

https://ja.wikipedia.org/wiki/%E7%89%B9%E7%95%B0%E5%80%A4%E5%88%86%E8%A7%A3

おすすめ記事

ガボールフィルター(Gabor Filter)

ガボールフィルターは、以下の式で表される。

gaussian(σ, z) * exp ( i * z )

ここで、gaussian項はガウス分布(正規分布)の関数、exp項のiは虚数を意味し、zは(x, y):位置情報 を表している。

上記のガウス分布は、2次元の場合

gaussian (σ, x, y ) = ( 1 / 4 * π * σ ) * exp ( -1 * (x2 + y2) / (4 * σ2) )

として表される。

上記ガウシアンの式をgnuplotで描画した例は以下の通り

σ = 1(標準偏差) として、gauss(x) = 1 / (2 * sqrt(pi)) * exp(- x * x / 4 ) をプロット(シグマを1としたため、式が簡略化)

ガボール gaussian(σ, z) * exp ( i * z )の分解

式には、exp(i * z)の虚数が入っている。このパートを、実部と虚部に分ける。(exp(ix) は cosx + i sinxと分解できる)

画像処理でガボールフィルターを用いて処理を施す場合は、この二つのパートに分離して作業をすることが多い。

ガボールフィルターには、ガウシアン項とExp項があり、Exp項には虚数があり、Cos項 と Sin項に分離できる。(exp(iθ) = cosθ + isinθ)(下記で説明する回転のsin, cosとは違うので混同しないように。)

それぞれの項目と、ガウス関数を掛け合わせたのは以下の通り。

ところで、本やWebでは、ガボールカーネルの関数は、上記で示した gaussian(σ, z) * exp ( i * z ) ではなく、

gaussian(σ, z) * (exp ( i * z ) – exp ( σ2 ) )

と表される。このexp ( σ2 )を引いている理由は、

gaussian(σ, z) * exp ( i * z ) を積分した結果が0にならない(直流成分が0とならない)ためである。(これはCos項:実数部が原因している。)

exp ( σ2 )を引いた画像は以下の通り。

Cos項とガウシアン関数の畳み込みカーネル:実部実数部カーネル(cos * gaussian) 3 – dimention

Sin項とガウシアン関数を畳み込みカーネル:虚部虚数部カーネル(sin * gaussian) 3 – dimention
splot cos(x) * gauss(x,y)
splot sin(x) * gauss(x,y)
gauss(x,y) = (1 / (4 * pi)) * exp(-1 * (x ** 2 + y ** 2) / (4))
splot cos(x) * gauss(x,y)
splot sin(x) * gauss(x,y)

として描画。(ガウス関数はσを1とおいた。そのため式が簡略化されている。)注意:実際に、直流をゼロにするためには、Cos項は、

      gauss(x,y) * ( cos(x) –  exp(-1) )

とすべき。

尚、カーネルにスケーリング(a倍)を施したい場合は、

 1 / sqrt(a) * gaussian( σ, z / a ) * (exp ( i * z / a ) – exp ( σ2 ) ) 

一般的に書くと、aj 倍したい場合、

 1 / pow(a, 1 / j) * gaussian( σ, z / aj ) * (exp ( i * z / aj ) – exp ( σ2 ) ) 

となる。

なお、回転を施したい場合は、zに回転行列をzにかければよい。

z = (x, y)とすると、

  [x’]     [ cosθ    sinθ ]  [x]

  [y’]     [ -sinθ   cosθ ]  [y]

z’ = (x’, y’)

参考文献

https://ja.wikipedia.org/wiki/%E3%82%AC%E3%83%9C%E3%83%BC%E3%83%AB%E3%83%95%E3%82%A3%E3%83%AB%E3%82%BF

おすすめ記事

正準相関分析(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

おすすめの記事