雑記ブログ【HELLO REALWORLD】始めました。

【PyTorch】多項式回帰MAP推定・変分ベイズ推論を試す

CAE・有限要素法のお勉強：ガラーキン法

KL reverseとKL forwardのメモ

変分推論とVAEまでの流れ

AutoMLとは？

変分モデルの書き方 Pyro

torch.jit を使ってみたのでメモ

Graph Neural Network 系のライブラリメモ

TensorFlow 2.0 のコードの書き方基本集（PyTorchとの比較）

はじめに

最近KerasからPyTorchに流れていく人たちが多く見受けられます。その中でて「Kerasで書いていたコードをPyTorchのコードで実装するにはどうすれば良いんだろう？」という声があります。要は、今まで使っていたフレームワークでやろうとしていたことを、別のフレームワークでやろうとしたときに、どんな対応になるのかが明確にわかっていればコードは圧倒的に書きやすくなるので、それを知りたいという話だと思われます。

ご要望に従って（？）、今後PyTorchからTF2.0に移りにあたり「PyTorchで書いていたコードに対応する、TF2.0でのコードはどんな感じになるんだ？」というところにお答えします。

線形回帰と学習のコード

データセット

下記のデータをPyTorchとTFで共通して使うことにします。

import numpy as np

# Hyper-parameters
input_size = 1
output_size = 1
num_epochs = 60
learning_rate = 0.001

# Toy dataset
x_train = np.array([[3.3], [4.4], [5.5], [6.71], [6.93], [4.168], 
                    [9.779], [6.182], [7.59], [2.167], [7.042], 
                    [10.791], [5.313], [7.997], [3.1]], dtype=np.float32)

y_train = np.array([[1.7], [2.76], [2.09], [3.19], [1.694], [1.573], 
                    [3.366], [2.596], [2.53], [1.221], [2.827], 
                    [3.465], [1.65], [2.904], [1.3]], dtype=np.float32)

PyTorch

import torch
import torch.nn as nn
import numpy as np



# Linear regression model
model = nn.Linear(input_size, output_size)

# Loss and optimizer
criterion = nn.MSELoss()
optimizer = torch.optim.SGD(model.parameters(), lr=learning_rate)  

# Train the model
for epoch in range(num_epochs):
    # Convert numpy arrays to torch tensors
    inputs = torch.from_numpy(x_train)
    targets = torch.from_numpy(y_train)

    # Forward pass
    outputs = model(inputs)
    loss = criterion(outputs, targets)
    
    # Backward and optimize
    optimizer.zero_grad()
    loss.backward()
    optimizer.step()
   

TF2.0

import tensorflow as tf
import numpy as np

L = tf.keras.layers

# Linear regression model
model = L.Dense(output_size)

# loss function
def loss_fn(y_predict, y):
    return tf.keras.losses.mean_squared_error(y_predict, y)

# optimizer
optimizer = tf.keras.optimizers.SGD(learning_rate=learning_rate)

for epoch in range(num_epochs):

    inputs = tf.convert_to_tensor(x_train)
    targets = tf.convert_to_tensor(y_train)

    with tf.GradientTape() as tape:
        y_predict = model(inputs)
        loss = loss_fn(y_predict, targets)
    grads = tape.gradient(loss, model.variables)
    
    # update prameters using grads
    optimizer.apply_gradients(zip(grads, model.variables))

違い

些細な違い：層の定義の仕方

PyTorchでは nn.Linear のコンストラクタにインプットの次元とアウトプットの次元を両方を渡してやる必要があります。一方でTFではL.Denseに対してアウトプットの次元のみを渡せばよく、インプットの次元は、最初にデータが入力されたときにそのデータの次元で自動的に決定されます。

この入力データによるインプットの次元を決定する恩恵は畳み込み層やLSTMなどのアウトプットを線形層にインプットするようなケースで感じるでしょう。

些細な違い：ロス関数の書き方

ここではロス関数の定義の仕方に違いがあります。PyTorchは色々なロス関数がnn.MSELoss()のようにクラスとして定義されているので、そいつのインスタンスを拾ってくるということになります。当然インスタンスの__call__メソッドに、具体的な計算が定義されているのですが、PyTorchは多くの場面で具体的な処理はメソッドとして隠蔽しつつ、内部でいろんな関数を呼び出すようになっています。

TF2.0ではtf.keras.lossesモジュールの関数を呼び出すという形式を取っています。今回、Python関数としてラッピングする意味は全く無いので、学習コードのところに直接 loss = tf.keras.losses.mean_squared_error(y_predict, y) と書いてしまってもいいでしょう。 基本的には適宜自分で必要な関数を呼び出すという書き方になります。

大きな違い：勾配計算とパラメータ更新

学習コードでその違いは顕著になっており、まずPyTorchではoptimizerに予め訓練するパラメータを渡してしまいます。そうすることで、optimizerは内部変数として訓練パラメータを保持するようになり更新時には更新メソッドであるupdate()を呼ぶだけで更新が可能になっています。

一方でTFではoptimizerは計算式を持っているだけで、どの変数をどの勾配を使って更新するのかは、後で更新メソッドapply_gradientsに対して引数として（訓練パラメータのリストとその勾配のリストをzipして）渡します。

勾配の計算にも違いが出ており、PyTorchは基本的に常に計算グラフを保持しながら計算を進めるので、lossは自分がどのように計算されたかを知ることができ、loss.backward()のようにbackward()メソッドを使うことで、lossを計算するにあたって用いられてtorch.tensorによる勾配を
各々のtorch.tensor の内部変数に格納
します。今回の場合は線形変換に使われるパラメータのtorch.tensorの勾配がそれぞれのパラメータの内部変数に入ります。前述の通り、optimizer は訓練するパラメータを保持しているので、それらの勾配にもアクセスができ、メソッドに引数を渡すことなく楽々と勾配計算から更新までができてしまうのです。

TFではwith tf.GradientTape() as tape コンテキストで囲んだ部分で計算グラフの履歴を残す計算をすることになります。このとき履歴は tape が保持するので、勾配計算などをする場合は tape に処理をしてもらうことになります。このコンテキスト内には色々な計算式が含まれうるので、tape.gradient()メソッドに微分される変数（今回はloss）と、微分する変数（model.variables）を明示的に渡してやります。逆にいろいろな計算が含まれていても、このように必要最低限の計算しか実施されません。

勾配を取り出したら前述の通り、optimizer の更新メソッド apply_gradientsに訓練パラメータと勾配の zipを渡してあげます。基本的に何を渡してどんな処理を行うのかをコードで明示して書いていくスタイルになります。

ニューラルネットワークの簡単な書き方

層をひたすら渡すだけで、簡単に多層のニューラルネットワークを書くことができるよう、それぞれSequentialというクラスが用意されています（おそらくkerasが最初に提供？）。

PyTorch

model = nn.Sequential(
    nn.Linear(in_size, hidden_size),
    nn.ReLU(),
    nn.Linear(hidden_size, hidden_size),
    nn.ReLU(),
    nn.Linear(1),
)

nn.Sequential のコンストラクタに層を順番に渡してあげます。

TF2.0

model = tf.keras.Sequential([
    L.Dense(hidden_size, activation="relu"),
    L.Dense(hidden_size, activation="relu"),
    L.Dense(1),
])

tf.keras.Sequential のコンストラクタに
層を順番に入れたリスト
を渡してあげます。

違い

見ての通り、似た書き方になりますが、渡すものが「複数の層」（PyTorch）なのか、「複数の層を格納した１つのリスト」（TF）なのかという違いがあり、地味ですがハマったりする人がたまにいます（これはKerasからPyTorchに行った場合にハマる可能性が在る部分かもしれません）。

また、TFでは層そのものが活性化関数を保持できるので行数を削減できます。ただし、バッチ正規化などを間に挟んだりドロップアウトを挟んだり、そういうケースにも対応できるように、層には活性化関数を持たせずに、活性化関数だけの層を挟む方法も用意されています。

model = tf.keras.Sequential([
    L.Dense(hidden_size),
    L.ReLU(),
    L.Dense(hidden_size),
    L.ReLU(),
    L.Dense(1),
])

これはKerasがTFにマージされ、PyTorchを意識し始めた段階で準備されたように思います。

畳み込みニューラルネットワーク

PyTorch

class ConvNet(nn.Module):
    def __init__(self, num_classes=10):
        super(ConvNet, self).__init__()
        self.layer1 = nn.Sequential(
            nn.Conv2d(1, 16, kernel_size=5, stride=1, padding=2),
            nn.BatchNorm2d(16),
            nn.ReLU(),
            nn.MaxPool2d(kernel_size=2, stride=2))
        self.layer2 = nn.Sequential(
            nn.Conv2d(16, 32, kernel_size=5, stride=1, padding=2),
            nn.BatchNorm2d(32),
            nn.ReLU(),
            nn.MaxPool2d(kernel_size=2, stride=2))
        self.fc = nn.Linear(7*7*32, num_classes)
        
    def forward(self, x):
        out = self.layer1(x)
        out = self.layer2(out)
        out = out.reshape(out.size(0), -1)
        out = self.fc(out)
        return out

TF2.0

class Cifar10Model(tf.keras.Model):
    def __init__(self):
        super(Cifar10Model, self).__init__(name='cifar_cnn')
        
        self.conv_block1 = tf.keras.Sequential([
            L.Conv2D(
                8, 
                5,
                padding='same',
                activation=tf.nn.relu,
                kernel_regularizer=tf.keras.regularizers.l2(l=0.001),
            ),
            L.MaxPooling2D(
                (3, 3), 
                (2, 2), 
                padding='same'
            ),
        ])

        self.conv_block2 = tf.keras.Sequential([
            L.Conv2D(
                16, 
                5,
                padding='same',
                activation=tf.nn.relu,
                kernel_regularizer=tf.keras.regularizers.l2(l=0.001),
            ),
            L.MaxPooling2D(
                (3, 3), 
                (2, 2), 
                padding='same',
            ),
        ])
        
        self.conv_block3 = tf.keras.Sequential([
            L.Conv2D(
                32, 
                5,
                padding='same',
                activation=tf.nn.relu,
                kernel_regularizer=tf.keras.regularizers.l2(l=0.001),
            ),
            L.MaxPooling2D(
                (3, 3), 
                (2, 2), 
                padding='same',
            ),
        ])
        
        self.flatten = L.Flatten()
        self.fc1 = L.Dense(
            256, 
            activation=tf.nn.relu,
            kernel_regularizer=tf.keras.regularizers.l2(l=0.001))
        self.dropout = L.Dropout(0.8)
        self.fc2 = L.Dense(10)
        self.softmax = L.Softmax()

    def call(self, x, training=False):
        x = self.conv_block1(x, training=training)
        x = self.conv_block2(x, training=training)
        x = self.conv_block3(x, training=training)
        x = self.flatten(x)
        x = self.dropout(self.fc1(x), training=training)
        x = self.fc2(x)
        return self.softmax(x)

違い

パディング

まず、畳み込み層の大きな違いはpadding の設定になると思われます。PyTorchではパディングを数値で渡しますが、TFの場合は'same'か'valid'で渡します。sameは縦横が変わらないようにパディングをするという意味で、validはパディングしないという意味になります。もしも縦横の大きさを変えたくない場合は same は非常に便利です。PyTorchの場合畳み込みのカーネルに応じて自分で計算することになります。

一方で、意図した大きさに縦横を明示的に変えていきたい場合は、TFの畳み込み層では実施することすらできません。代わりに畳み込み層の方ではvalidでパディングをしない設定にしつつ、L.ZeroPadding2Dという層を使うことで明示的にパディングを実施してやることができます。

畳み込み層→線形層

畳み込み層から線形層に変わる部分ではテンソルの形を調整する必要があります。 PyTorchでは自分で形を変えてあげる必要がありますが、TFには線形層に入れるためのL.Flatten()があるため、そのようなコードは省くことができます。ただし、結局 tensor.reshape(batchsize, -1) としてやればよく、batchsize=tensor.shape[0] などで簡単にその場で取得できるので大したコードの節約にはなりません。

このL.Flatten()が活躍するのは、実はtf.keras.Sequential を利用するようなケースです。

traininigフラグ

ここでついでに解説しますが、PyTorchのtorch.nn.Moduleにあるmodel.train()やmodel.eval()のようなメソッドはtf.keras.Modelにはありません。callメソッドで訓練時とテスト時の振る舞いを明示的に分ける必要があります。ちなみにtf.keras.Sequentialクラスはcallメソッドを自分で書かないのですが、ちゃんと最初からtraining引数を取れるようになっており、L.Dropoutなどに自動で適用されるようになります。

RNN

今回はBidirectional RNNをLSTMを使って見ていきます。

PyTorch

class BiRNN(nn.Module):
    def __init__(self, input_size, hidden_size, num_layers, num_classes):
        super(BiRNN, self).__init__()
        self.hidden_size = hidden_size
        self.num_layers = num_layers
        self.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True, bidirectional=True)
        self.fc = nn.Linear(hidden_size*2, num_classes)  # 2 for bidirection
    
    def forward(self, x):
        # Set initial states
        h0 = torch.zeros(self.num_layers*2, x.size(0), self.hidden_size).to(device) # 2 for bidirection 
        c0 = torch.zeros(self.num_layers*2, x.size(0), self.hidden_size).to(device)
        
        # Forward propagate LSTM
        out, _ = self.lstm(x, (h0, c0))  # out: tensor of shape (batch_size, seq_length, hidden_size*2)
        
        # Decode the hidden state of the last time step
        out = self.fc(out[:, -1, :])
        return out

TF2.0

class BiRNN(tf.keras.Model):
    def __init__(self, hidden_size=10, num_layers=2, num_classes=10):
        super(BiRNN, self).__init__(name='mnist_rnn')
        self.hidden_size = hidden_size
        self.num_layers = num_layers
        
        self.lstm = self.get_lstm_layers(hidden_size, num_layers)            
        self.fc = L.Dense(num_classes, activation="softmax")
    
    @staticmethod
    def get_lstm_layers(hidden_size, num_layers):
        lstm_layers = []
        # we need all sequence data. write return_sequences=True! 
        for i in range(num_layers-1):
            lstm_layers.append(
                L.Bidirectional(
                    L.LSTM(units=hidden_size, 
                                         return_sequences=True)
                )
            )
        # the final layer return only final sequence
        # if you need all sequences, you have to write return_sequences=True.
        lstm_layers.append(
            L.Bidirectional(
                L.LSTM(units=hidden_size)
            )
        )
        return tf.keras.Sequential(lstm_layers)
        
    def call(self, x):        
        # Forward propagate LSTM
        out = self.lstm(x)
        out = self.fc(out)
        return out

違い

大きな違い：多層化

まず、PyTorchのnn.LSTMは最初から多層にすることを想定しており、num_layersなる引数を持っています。TFで多層化したい場合は、SequentialにLSTMがたくさん入ったリストを渡してあげる必要が在るため、この辺で面倒なコードを書く必要が出ます。

些細な違い：Bidirectional

PyTorchのLSTMではBidirectionalにしたければコンストラクタでその指定が可能です。TFではL.Bidirectoinalでラッピングすることになります。

大きな違い：戻り値の並び

PyTorchのLSTMは (seq_len, batch, features) という形でテンソルが返ってきます（それプラス隠れ状態とセル状態も帰ってきます）。一方でTFでは(batch, seq_len, features) という形でテンソルが返ってきます。ただし、PyTorchは(batch, seq_len, features)で返ってくるように指定することも可能です。上記のコードではPyTorch側をTF側に揃えています。

もしも最後のシーケンスのアウトプットのみが欲しい場合は、PyTorchの場合はスライシングで対応することになりますが、TFの場合LSTMのコンストラクタに全てのシーケンスを返させるか、最後のシーケンスのみを変えさせるかを指定できます。スライシングもできるので、特に迷ったのであれば全て返させても良いかもしれません。ただ、一旦全てのシーケンスを何らかの変数に保存して、再度取り出すのは無駄と言えば無駄かもしれません（そう影響は無いだろうが）。

学習

学習に関しては既に線形回帰の方で比較をしたのですが、実はTFの方はモデルを作った後 model.compileしてmodel.fitをすればKerasの学習ラッパーをそのまま使えます。凄く便利です。

今回は疲れたのでこの辺で。

github.com

確率的プログラミング言語Pyroと変分ベイズ推論の基本

はじめに

タイトルの通り基本だけ書きます。 Pyroは解説できるほど触っていないので、大したことは書けませんが、何も知識が無いよりはとっつきやすくなるであろうことを書いておきます。

ベイズ推論

モデリング

ベイズモデリングとは観測データ $X$ を、パラメータだとか潜在変数だとか呼ばれる「未観測の確率変数」を使ってモデリングする試みです。例えば $X$ が正規分布に従うと思っているのならばパラメータ $\mu, \sigma$ を用いて、$ N(X | \mu, \sigma) $という確率モデルを考えてみたりします。

一般に $p(X | Z)$ とか書いておくことが多いです。 $Z$ というのはパラメータやら潜在変数やら、未観測の確率変数で$X$ の方が観測した手元のデータになります。ベイズモデリングでは未観測の確率変数 $Z$ に対して事前分布 $p(Z)$ を考え、確率モデルと事前分布のタプルで $(p(X|Z), p(Z))$をベイズモデルとして扱います。

事後分布

ベイズモデルを決めると、ベイズの定理から下記の事後分布を書き下すことが出来ます。

$$ p(Z | X) = \frac{p(X|Z)p(Z)}{p(X)} $$

分子は自身が仮定したベイズモデルによって構成されています（それが容易に計算できるかはさておいて）。分母は観測したデータ $X$ によってのみ構成されているので、すでに確定値です（確定値ではあるが計算できるかはこれも保留）。

このように未観測の確率変数 $Z$ に関する分布を観測済のデータ $X$ の条件付き分布（事後分布）として求めることをベイズ推論と言います。

予測分布

さて、 $Z$ とは具体的には仮定した正規分布の平均パラメータかもしれませんし、カテゴリ分布の確率を表すパラメータかもしれませんし、はたまた、物理的実体と照らし合わせた変数かもしれません。事後分布を求めることで、手元に在るデータから自分が仮定した未観測の確率変数に関する分布を得ることができるわけですが、私達が興味を示すのは、今後手に入る $X$ に関する予測であるケースが多いです。

あくまで $Z$ という未観測の確率変数は、 $X$ をモデリングするための仮定として置いたものであり、その潜在変数そのものよりも、置いた仮定によってどれだけ $X$ を上手く表せるかが重要になってくるケースが多いということです（というかそういう汎化性能の評価をすべきだと思う）。

ベイズで使われる予測分布は下記の式で表されます（すでに持っているデータ $X$ を使って $x_{new}$ の分布を表す）。

$$ p(x_{new} | X) = \int_Z p(x_{new}|Z)p(Z|X)dZ $$

これはすなわち、もともと仮定したベイズモデル $(p(X|Z), p(Z))$の確率モデル $p(X | Z)$ の $X$ の部分をこれから予測したい $x{new}$ に置き換え、 p(x{new} | Z) という確率モデルを事後分布 $p(Z|X)$ で重み付けして分布を構築するということです。

実際に使われる予測分布

多くのケースでは

$$ Z_i \sim p(Z|X) $$

とサンプリングを行って、

$$ p(x_{new} | Z_i) $$

という確率モデルを得る。ということをできるだけ大きな数 $N$ 回繰り返すことで、

$$ \frac{1}{N} \sum_{i=1}^N p(x_{new} | Z_i) $$

と $N$ 個の確率モデルの和で予測分布を近似して使います。 $Z_i \sim p(Z|X)$ は当然事後分布 $p(Z|X)$ により選ばれやすい $Z$ が頻出するはずですので、それに応じて選ばれやすい $Z_i$で構築される確率モデル $p(x_{new} | Z_i)$ の割合が多くなるので、元々の予測分布を近似できるという仕組みです（PPLでは数式として解くのが困難な積分はこの方法によって近似的に処理されることが頻繁にあるので、この部分を納得しているだけでだいぶ感覚が違うはずです）。

Pyroの基本

結局のPyroを始めとしたPPLでやろうとしていることは、

• データ $X$ に対して確率モデル $p(X|Z)$ を作る
• $Z$ に対して事前分布 $p(Z)$ を割り当てる
• 事後分布 $p(Z|X)$ を求める
• 予測分布 $p(x_{new} | X )$ を得る

という操作を手軽に実行する仕組みを与えたいというものです（もちろんもっとプリミティブに、確率変数というものを上手に扱う仕組み、例えば微分可能にする仕組みなどを入れることも重要な仕事のひとつである）。

それぞれの操作がどのようなAPIになっているのかを把握できれば、それだけでPPLの見通しがグッと良くなります。今回はPyroのAPIを見ていきましょう。

Pyroの確率変数の取扱

pyroでは確率変数は一貫してpyro.sample('name', distribution_fn)によって表すと決められています。では確率変数が観測データであるのか潜在変数であるのかをどのように見分けるのかというと、pyro.sample('name', distribution, obs=observed_data) として観測データを渡しておいたものだけが観測データとして扱われます。

単に distribution.sample() メソッドを利用した場合、それはPPLとしての確率変数という扱いはされず、言わば普通に numpy.random.randnなどを利用したのと同じになってしまいます（リッチな擬似乱数モジュールに過ぎないということ）。

pyro.sample()という関数は非常に重要なので必ず抑えましょう。

Pyroのハイパーパラメータの取扱

ハイパーパラメータのように固定値として使う物は普通にPythonの数値として、あるいはtorchにTensorとして与えれば問題ありません。個々に関しては普通のプログラミング同様なので困惑するところはないでしょう（困惑するとしたらベイズモデリング自体に困惑している可能性があります。ベイズモデリングでも決め打ちになるハイパーパラメータはあります。いくつか決め打ちを試してみて、その中で一番良さげなのを選ぶという手法もありますが）。

Pyroでの変分パラメータの取扱

Pyroで変分パラメータを扱う場合は、pyro.param("name", tensor)を用います。こいつで宣言された変数tensorは変分パラメータとしてPyroが認識するようになります。変分推論で事後分布を近似する際には、パラメトリックに置いた事後分布のパラメータを上記の関数で宣言することを忘れてはいけません。

変分ベイズ推論のコード：確率モデル

基本的なpyroのパーツとしては、本当に pyro.sample と pyro.param が分かれば一先ず…という感じです。あとは pyro.distributions を使うのですが、これは torch.distributions をラップしてる pyro用の確率分布モジュールで、特に特有の難しさはありません。

まず、ベイズ確率モデルは確率モデル $p(x|\theta)$ と事前分布 $p(\theta)$の組でできるのでした。ここで、事前分布 $p(\theta)$ は何らかの確率分布を仮定するのですがから、その確率分布を表すためのパラメータ $\alpha$ （俗にいうハイパーパラメータ）を決め打ちで準備する必要が在ることに注意してください。

すなわち、正確には

$$ p(x \mid \theta) $$

という確率モデルと

$$ p(\theta \mid \alpha) $$

という事前分布（$\alpha $ は決め打ち）を仮定することから始めます。

ここでは、データ $x$ に対する確率モデルとしてベルヌーイ分布 $p(x \mid f)$ と、事前分布としてベータ分布 $p(\f \mid \alpha, \beta)$ を選ぶこととしましょう。ハイパーパラメータは事前分布のパラメータ数に応じて数が変わることと、これに関しては単なる決め打ちの値になることに注意して下記のコードを読みましょう。

def model(data):
    alpha0 = torch.tensor(10.0)
    beta0 = torch.tensor(10.0)
    # sample f from the beta prior
    f = pyro.sample("latent_fairness", dist.Beta(alpha0, beta0))
    # loop over the observed data
    for i in range(len(data)):
        # observe datapoint i using the bernoulli
        # likelihood Bernoulli(f)
        pyro.sample("obs_{}".format(i), dist.Bernoulli(f), obs=data[i])

上記のコードを数式で表すならば

$$ \alpha = 10 $$ $$ \beta = 10 $$ $$ f \sim \rm beta(f \mid \alpha=10, \beta = 10) $$ として、$f$ を使って下記のサンプリングを実施。 $$ x_i \sim {\rm Bern}(x\mid f) $$

ということになります。これによって

$$ {\rm Beta}(f\mid \alpha, \beta) \prod_{i}^{N} {\rm Bern}(x \mid f) $$

のような計算が実施されたことになるのです（実際には内部的にはこの対数を取って、対数同時確率を計算しているであろう）。ちなみに観測データである data は $model$関数の引数として渡して、内部の対応する確率変数を pyro.sample関数で生起させるときにobs 引数で紐づけておくことになります。

この辺りがpyroを扱う際の肝になってきます。 ちなみに今回はわかりやすく forループを用いていますが、 with pyro.plate が処理の並列化まで自動で行ってくれるため高速化に使えます（Stanでいうところのベクトル化）。

変分モデル

上記の確率モデル（正規化はされていないが）

$$ {\rm Beta}(f\mid \alpha, \beta) \prod_{i}^{N} {\rm Bern}(x \mid f) $$

に対して、事後分布 $p(f \mid x)$ が我々の知りたいものになります。今回、実はベルヌーイ分布に対してベータ分布を事前分布として選択しているため、解析的に事後分布を求めることができますが、今回は練習のため変分推論を行ってみます。

変分推論のポイントは、事後分布が解析的に解けなかったり、計算が面倒であったり、サンプリングすら計算量的によろしくなさそうな場合に、事後分布を適当なパラメトリックな関数（確率分布の形状として考えられそうなものを選ぶ）として置いてしまうというところです。

$$ q(f) = {\rm Beta} (\alpha', \beta') $$

というパラメトリックな関数として置いてみましょう。もちろんベータ関数で置かなくても良いです。事後分布は正規分布っぽいんじゃないか！？と思うならそれでも良いですし、実は減衰する正弦波みたいな形なんじゃないか！？と思うならそうしてください。そのような形状の者の中で、ELBOが最大化されるような形状が選ばれます。

今回は、解析的に解くとベータ関数で表されることを知っているのでインチキをして事後分布をベータ関数としておいてしまいます。

def guide(data):
    # register the two variational parameters with Pyro.
    alpha_q = pyro.param("alpha_q", torch.tensor(15.0),
                         constraint=constraints.positive)
    beta_q = pyro.param("beta_q", torch.tensor(15.0),
                        constraint=constraints.positive)
    # sample latent_fairness from the distribution Beta(alpha_q, beta_q)
    pyro.sample("latent_fairness", dist.Beta(alpha_q, beta_q))

ここで、いま置いたベータ関数（で表されるベータ分布）によって $f \sim \rm Beta$ と出てくると決めたのであれば、それで表現される形の中でそれらしい事後分布が選ばれるだけで、それ以外の物を見つけることはできなくなります。今回は、事後分布の形状の仮定としてベータ分布を置いたことは（正解を知っているので）非常に筋の良い選択ということになりますが、通常はこの部分は手探り近いことを認識しておきましょう。

また、alpha_q と beta_q は変分パラメータなのでpyro.param で変数を作ることになります。引数は名前と初期値と、制約です（今回は正の値という制約が入っている。これは分散などをパラメトライズする場合にもしばしば使うので覚えておくこと）。

pyroの注意点としては、model が確率モデルと事前分布を記述し、 guide が変分モデルを記述すると明確に分けつつも、guide の引数をmodel と合わせる風習が在るようです。当然今回の場合 $q(f) $という変分モデルはデータを引数に取っていない単なるパラメトリックな関数としてなので、内部でデータは使われていないことに注意しましょう（ところがどっこい、VAEなどの場合、ニューラルネットワークのパラメータ自体が全て変分パラメータで、 decoderを確率モデル $p(x|z)$ 正規乱数を事前分布 $p(z)$ encoderを変分モデル $p(z | x)$ として推論を行うので、guide側もdata を内部で使う。 この辺はVAE自体を理解していない場合に混乱をするもととなるので注意が必要。そもそもVAEはパラメータを単にKL損失関数とした普通の勾配学習と思った方が実装は楽である。）。

学習コード

学習は驚くほど簡単です。

# set up the optimizer
adam_params = {"lr": 0.0005, "betas": (0.90, 0.999)}
optimizer = pyro.optim.Adam(adam_params)

# setup the inference algorithm
svi = pyro.infer.SVI(model, guide, optimizer, loss=pyro.infer.Trace_ELBO())

n_steps = 5000
# do gradient steps
for step in range(n_steps):
    svi.step(data)

と記述するだけ！びっくりだね！！ SVIなるクラスが確率モデル model と変分モデル guide を渡せば学習をやってくれるというスグレモノになっています。 損失関数はELBOが使えますが、今後他のものも追加されていくのかどうかは謎…。

とりあえずモデルの設計者側が実施するのは基本的に、確率モデルを書いて、変分モデルを書いて、あとは変分推論用のクラスにそいつらを渡してやるという作業だけです。簡単ですね。もちろん確率モデルを階層ベイズモデルにしたければそうすればいいでしょう。

変分推論のカスタマイズ

変分推論のクラスを使うと学習がまるっきり隠蔽されます。少し低レベル側も見ましょう。

# define optimizer and loss function
optimizer = torch.optim.Adam(my_parameters, {"lr": 0.001, "betas": (0.90, 0.999)})
loss_fn = pyro.infer.Trace_ELBO.differentiable_loss
# compute loss
loss = loss_fn(model, guide)
loss.backward()
# take a step and zero the parameter gradients
optimizer.step()
optimizer.zero_grad()

というコードがピンと来る方は、まさしくpyroがpytorchをラッピングしたものであるとすぐに分かるはずです。 変分推論は変分パラメータを普通の学習パラメータと思った場合に、勾配学習で最適化を行っているに過ぎません。ただし、その時の損失関数が確率モデル（＋事前分布）と勝手に決めた計算しやすい変分モデルとのKLダイバージェンスないしELBOで規定されているという代物になっています。

もちろん学習後は、変分モデル（事後分布の近似関数）を使って、ベイズ予測分布を構築したりすることができる点で、単なる確率モデルのパラメータの点推定とは異なるのですが、事後分布をパラメトライズして最適化を行うという形式は、計算テクニック上同じになってきます（ここらへんで混乱する人は、まずベイズ推論の基本を勉強する必要があるような気がします。かくいう自分もベイズに触れ始めたときは大層混乱した）。

更に損失関数として使っている pyro.infer.TraceELBO.differentiable_loss もよくよく見てみると

# note that simple_elbo takes a model, a guide, and their respective arguments as inputs
def simple_elbo(model, guide, *args, **kwargs):
    # run the guide and trace its execution
    guide_trace = poutine.trace(guide).get_trace(*args, **kwargs)
    # run the model and replay it against the samples from the guide
    model_trace = poutine.trace(
        poutine.replay(model, trace=guide_trace)).get_trace(*args, **kwargs)
    # construct the elbo loss function
    return -1*(model_trace.log_prob_sum() - guide_trace.log_prob_sum())

svi = SVI(model, guide, optim, loss=simple_elbo)

こんな感じで（当然と言えば当然だが）数式レベルで確率モデルと変分モデルの対数確率の和（の負）を計算してることが丸見えになります。計算に使われている poutineモジュールやその他諸々のメソッドについては公式Docに任せますが、兎にも角にもpyro内でインスタンスを管理する上で色々上手に隠蔽しているということです。

pyroについて

仮に数式レベルで書くことに抵抗がそれほど無いのであれば、TensorFlow Probabilityをおすすめします。 個人的には pyro はTensorFlowで言うところのKeras的小回りの効かなさを感じます。やりたいことがシンプルであればpyroの枠組みで済みますが、損失関数レベルでごちゃごちゃ触るなら、隠蔽されている中身を紐解くのが帰って面倒です。

統計的学習理論でのコトバの定義が分かりやすかったのでまとめ

はじめに

青いシリーズの統計的学習理論

統計的学習理論 (機械学習プロフェッショナルシリーズ)

こいつに書かれているコトバの定義がわかり易かったので、まとめておきます。 機械学習が何たるものであるかを端的に表現することができるので、共通の認識になると良いな。

コトバの定義

データ

観測によって得られる情報。既知の情報だけでなく将来に渡って得られる情報のことも指します。 一般には数値で表現されることが多く、多くの場合はベクトルや行列で表現されます。近年はもっと多次元の配列（テンソルと広く呼ばれている）としてデータを表すようになってきました。

学習データ

データの中でも機械学習モデルの学習に利用するデータを「学習データ」と言います。あるいは「訓練データ」とか「training data」などとも呼ばれます。

検証データ

学習途中段階、あるいは学習後のモデルの性能評価に使われるデータです。「validation data」と表現される場合もあります。多くの場合は機械学習モデルの正則化パラメータ等、学習の際に人間が設定しなければならないハイパーパラメータと呼ばれる数値の選定に使われます。

テストデータ

学習後の機械学習モデルの汎化性能を見積もる際に使われます。想定している状況は、「将来手に入るであろうデータに対しても上手く対応できるモデルであるか否かを試す」という立場であるため（もちろん、試す以上は既に手に入っているデータを使うのだが…）、予め学習データと検証データの何れからも完全に除外されたデータを使う必要があります。

テストデータの結果を真髄に受け止めず、このテストデータに対してモデルの調整をし始めてしまったら、テストデータに対して過剰適合したモデルが出来上がってしまうので注意が必要です（それで良いモデルができたと思っても糠喜びである）。

入力データと出力データ

機械学習モデルに入力するデータのことを入力データと呼び、機械学習モデルから出力されるデータのことを出力データと呼びます。入力データ $x$ で出力データ $y$ のときデータをこれらの組で $(x, y)$ と表します。特に $y$ が離散的で有限集合に値を取る場合はラベルと呼ばれます。 $y \in \{0, 1\}$ というケースでは特に2値ラベルと呼ばれ、要素の数が 3以上になると多値ラベルと呼ばれます。

入力 $x$ が取りうる空間を入力空間 $\mathcal X $ と呼びます。出力 $y$ が取りうる集合を $\mathcal Y$ と表します。

機械学習モデル

仮説

入力空間 $\mathcal X$ から出力集合 $\mathcal Y$ への関数

$$ f : \mathcal X \rightarrow \mathcal Y $$

を仮説と呼びます。急に抽象的になりましたが、要は何らかの機械学習モデルのことを指しています。下記の仮説集合と加えて考えると非常に分かりやすいはずです。

仮説集合

仮説 $f$ の集合を仮説集合 $\mathcal F$ と呼び、機械学習のアルゴリズムとは仮説集合 $\mathcal F$ の中から適切な仮説 $f$ を選び出す処理と言えます。

例：入力 $x \in R^d$ の機械学習モデルを $f(W, x) = Wx $ と置いた場合、 仮説集合は行列 $W$ で表現できるような線形変換全体を表し、

$$ \mathcal F = \{x \mapsto Wx \mid W\in R^{D \times d} \} $$

という集合の中から何らかの学習アルゴリズムによって具体的に $W$ を定めるということになります。

ニューラルネットワークの構造を決めたり、SVMでカーネルを定めたりすることが仮説集合を決めること（どのような仮説が含まれるか、学習結果としてどんなものが選ばれうるかを規定）に相当し、学習アルゴリズムを走らせることでパラメータを定め、学習結果を得ることが仮説集合の中から仮説を一つ選び出す操作になります。

「仮説」という表現はなかなか良いなと思っていて、確かに学習結果は正しい姿を保証しているわけでも何でも無い、有限データと勝手な仮定のもとで導いただけのものであるので、それを端的に表しているなと感じています。

損失

損失関数

損失関数とは、端的に言えば仮説の良し悪しを測るための尺度です。特に値が小さいほど、良い仮説であると言えるような尺度を導入します。例えば、 仮説を $f : \mathcal X \rightarrow \mathcal Y$ が入力 $x$ に対して $y'_i = f(x_i)$ という出力をするとしましょう。実際の入出力データが $(x_i, y_i)$ であるとすれば、下記のような損失関数を選ぶことができます。

$$ \rm {loss} (y', y ) = \sum_i (y'_i - y_i)^2 $$

という計算式は実際の出力データと仮説の予測出力との二乗誤差を損失関数と定めたことになります。

損失関数はどうすれば良い仮説を選んだと言えるのかを考え、機械学習モデルの設計者が選ぶものであることに注意してください。また、計算の都合上で本来測りたいはずの尺度を取り扱うのが困難な場合（例えば微分ができないなど）、都合の良いような関数に少しだけ変更しつつ本来の目的が失われないような改造を加えたりします。そういう場合には「代理損失」などと強調したコトバが使われるケースもあるようです。

予測損失

機械学習の目的は適切な損失関数のもと、未知の新たなデータ（テストデータ）に対する損失が小さくなるような仮説を手に入れることです。 ここで、未知の新たなデータに対する損失を「予測損失」と言います。ここで仮説 $f$ に対する、予測損失 $R(f)$ はテストデータ $(X, Y)$ の分布を $p(X, Y)$ として、下記で定義されます。

$$ R(f) = \mathbb E_{p(X, Y)} [{\rm loss}(f(X),Y)] $$

すなわち、テストデータの損失のテストデータの分布に関する期待値になります。テストデータは未知の新たなデータを想定しており、当然、我々はその分布に関する完全な情報を持ち合わせていないことに注意しましょう。

経験損失

一方で、訓練データが手元に $N$ 個あり、各１つ１つのデータ $(X_i , Y_i)$が $1/N$ の確率で生起するような（実に都合の良い）ケースを考えると、そのような訓練データの分布 $q(X, Y)$ として、

$$ \hat R(f) = \mathbb E_{q(X, Y)} [{\rm loss}(f(X), Y)] = \frac{1}{N} \sum_i {\rm loss}(f(X_i), Y_i) $$

と表すことができます。これを経験損失と言い、通常は経験損失を用いて学習を行います。経験損失と予測損失と最たる違いは、期待値を計算する際の分布であり、もしもテストデータと訓練データの分布が全く同じものであるとするならば、経験損失と予測損失は期待値の意味で一致し、不偏推定量となります。

変分ベイズ法の心２

はじめに

下記記事の続きで、お気持ちは理解している前提で進みます。

www.hellocybernetics.tech

変分ベイズ法の戦略

基本の復習

データ $D = {x_1, \cdots, x_N}$ が手元にあるときに確率モデル $p(x|\theta)$ と事前分布 $p(\theta)$ を設計するのがモデリングの第一歩でした。するとベイズの定理（あるいは乗法定理）から、下記の事後分布を獲得することがベイズ推論の目標になります。

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

さて、この分布を推論したときに最も嬉しい結果は $\theta$ を表す真の分布 $\hat p(\theta) $ にピッタリ一致することでしょう。しかし、私達は真の分布を知らないので、どのようにこれに近づけて行くかというのが問題になります。

一つの解決法として、真の分布に収束するような確率遷移核を上手に設計し、その確率遷移核を持つマルコフ連鎖でサンプリングを繰り返し、ヒストグラムを描かせることで分布を得るMCMCがあります。

もう1つが今回紹介する変分ベイズ法です。変分ベイズ法では推論したい分布をパラメトリックな関数として仮定し、適切な評価指標のもと仮定した近似分布の出来栄えを評価しながら、パラメータを少しずつ変えていく手法です（推論が最適化問題に帰着されている）。

すなわち、事後分布が

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

と表されており、こいつの分布の形状を表す関数を $q(\eta)$ と書き表すことにします。$\eta$ を上手く調整して事後分布を近似してやりたいという話です。実際にはこいつは $\theta$ を発生させる分布であるので、 $q(\theta ;\eta)$ と書いたり、 $q(\theta | \eta)$ と書いたり $q(\theta)_{\eta}$ と書いたりする場合もあります。この $q$ は兎にも角にも $\theta$ を発生させる分布であり、私達がこれから調整するのは $\eta$ であるということを明確に意識しておきましょう。

分布の評価指標

近似分布の出来栄えを評価するのに使われるのが、KLダイバージェンスです。天下り的ですが、KLダイバージェンスは確率分布 $q$ と 確率分布 $p$ の離れ具合を評価することが出来ます。そこで下記を最小化することで、 $q(\theta; \eta)$ の $\eta$ をいろいろ調整してみて $p(\theta | D)$ に最も近い分布を見つけてやろうと試みます。

$$ KL [q(\theta; \eta) : p(\theta | D)] = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta $$

今、どういう状況になったのかを振り返ります。確率分布 $p(\theta |D )$ を求めたいというのがベイズ推論でした。そこで、こいつを一先ずでたらめな $q (\theta)$ で表して少しずつ分布の形を変えながら良さげなものを探そうと考えます。ところがこれでは抽象的すぎるし、分布の変え方をどうするのかもよくわからんので $q(\theta ; \eta)$ とパラメトリックな関数にして、 $\eta$ の調整に問題をすり替えました。そして、どんな分布の形が良いのか（どんな $\eta$ が良いのかを）をKLダイバージェンスで図ろうと決めたのです。

KLダイバージェンスは見た目が複雑に見えますが、もはや $\eta$ の関数に過ぎません。

ただし、ここで注意が必要です。私達は今 $q(\theta ; \eta)$ と事後分布の形状を制限してしまいました。パラメトリックに表された関数で表現できる形状以外のものは $\eta$ をどれだけ調整したところで獲得できないことになります。しかし、それでも、表しうる形状の中ではKLダイバージェンスの意味で最も近い分布を得ることが出来ますので一先ず良いとしましょう。

もう1つ注意すべきは、パラメトリックに関数を置いたことで最適化の進み方も決めてしまっている点です。何度も申してきたとおり、分布の形状を知りたいというのが問題であってその形状を調整するために $\eta$ というパラメータを導入したのでした。言わばこれは分布のイジりやすさを向上するために導入したものであって、それ以外の意味は基本的にはありません（仮定した分布によほどの意図がない限り、こいつはただの調整役に過ぎない）。

もしもパラメータの置き方をもっと工夫してみたり、あるいはパラメータを置かずとも近似分布 $q$ を調整する手段があるならばそれでも構わないのです。よほど良い $q$ の選び方をすれば、上手に調整する方法が何かしら見つかるかもしれませんね。

ELBO

引き続いて、KLダイバージェンスを最小化することをもう少し詳しく見ます。

$$ \begin{align} KL [q(\theta; \eta) : p(\theta | D)] & = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta \\ & = \int q(\theta; \eta) \log \frac{q(\theta; \eta) p(D)}{p(\theta, D)} d\theta \\ & = \int q(\theta; \eta) \log p(D)d\theta + \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta, D)} d\theta \\ & = \int q(\theta; \eta) \log p(D)d\theta - \int q(\theta; \eta) \log \frac{p(\theta, D)}{q(\theta; \eta)} d\theta \\ \end{align} $$

ひたすら式変形を実施しました。ここで第一項に関しては $q(\theta; \eta)$ は確率分布であるので $\log p(D)$ と表すことができます。これは対数周辺尤度、あるいはエビデンスと呼ばれる値です。引数がすでに手元にある確率変数の実現値 $D$ であるので、$p(D)$ は具体的に値を計算できる定数にすぎません（ただし $p$ の形を私達は知らないので今は保留…）。

第二項の方がEvidence Lower Bound、略してELBOと呼ばれるものになります。ELBOは名前の通り、エビデンス（対数周辺尤度）の下界を示しており必ず、

$$ \log p(D) \geq \log \frac{p(\theta, D)}{q(\theta; \eta)} d\theta $$

となることが知られています（話は逆で、そうなることがわかったのでELBOと呼ばれている）。このエビデンスとELBOが答式で結ばれる条件が、

$$ KL [q(\theta; \eta) : p(\theta | D)] = 0 $$

となること、すなわちベイズ推論の結果、$q$ が 事後分布 $p$ をピタリと当てられた時になります。なので変分推論を行う際には「ELBOを最大化する」という表現が使われることも多々あります。エビデンスの方は直接は計算できない保留中の定数でしたので、 ${\rm KL} = {\rm evidence} - {\rm ELBO}$ となっており、KLを最小化することと ELBOを最大化することは同じことをしていることになります。

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

$p(D)$ はデータと $p$ が決まれば自動的に定数となります。当然データは手元にあるものを使うのでこの時点で確率変数の確定値となっており変数として扱う必要はありません。言わば単なる正規化定数（ちゃんと確率の和が1になるように割り算するための定数）のように振る舞います。事後分布の形状を決めているのは分子の方です。

分子の形状を上手く近似してやることができれば求めたい事後分布が求まりそうです。その分子の形状に相当するものを $q(\eta)$ という関数で近似することにします。

$$ q(\eta) = p(\theta | D) = \frac{p(D|\theta) $$

変分ベイズ法の具体的手段

関数 $q$ をどのように置くのか

さて、上記までの話で兎にも角にも

$$ KL [q(\theta; \eta) : p(\theta | D)] = \int q(\theta; \eta) \log \frac{q(\theta; \eta)}{p(\theta | D)} d\theta $$

を最小化しましょうという話になりました。 さて $q(\theta ; \eta)$ はどのように置くと良いのでしょうか。一つは正規分布のような分かりやすい分布を置いて、 $\eta = (\mu, \sigma)$ みたいにしてしまうことです。そうすると、実際の $\theta$ の事後分布が正規分布でなかったとしても、正規分布の中で最も実際のサンプルに近いものを探すことになります。

もう１つは

$$ \q(\theta, \eta) = q(\theta_1, \eta_1)\cdots q(\theta_K, \eta_K) $$

と置いてみることです（ここで求めたい $\theta$ が$K$次元ベクトル、あるいは $K$個のパラメータを持っていたとしています）。これは各々のパラメータが（データの条件付き）独立であることを仮定していることになります。独立性を仮定しているので、各パラメータ間の相関はつかめなくなることに注意してください。仮に $\theta$ が行列ならば、 各成分それぞれを独立とみなしても良いですし、それがやりすぎだと思ったら、各列が独立とみなして各列が多次元正規分布であると仮定するような方法もあるでしょう。

このように複数のパラメータを有する分布を、個々のパラメータの単純な分布の積に分解してしまう方法を平均場近似と言います。

平均場近似は、事後分布を推論する際に「うまい分布の変え方」を上手くひねり出すときに使えます。例えば、$\eta$ を変分パラメータとして分布の調整役を担ってもらう方法を集中的に取りざたしてきましたが、単純な分布の集まりに変えてしまえば、うまい分布の調整の仕方をパラメータに頼らず考案できる可能性もあるのです（分布をイジる主たる手法が $\eta$ の導入であるというだけで、実は必ずしも使わなくていい…）。

個々らへんをゴリゴリに手計算しながら学びたい場合は

機械学習スタートアップシリーズ ベイズ推論による機械学習入門 (KS情報科学専門書)

をおすすめします。

変分ベイズ法の心

ベイズ推論の基本

ベイズモデリングの概要については下記の記事を参考にしてください。

www.hellocybernetics.tech

概要をさらっとなぞると、ベイズ推論の基本的な話としては、観測データ $x$ の真の確率分布 $\hat p(x)$ を知る由もないので、確率モデル $p(x | \theta)$ でモデル化し、更にパラメータ $\theta$ にも事前分布 $p(\theta)$ を仮定します。

$$ p(x, \theta) = p(x | \theta)p(\theta) = p(\theta | x) p(x) $$

という確率分布に対していつでも成り立っている乗法定理から、

$$ p(\theta | x) = \frac{p(x|\theta)p(\theta)}{p(x)} $$

とできます。そこで 各 $i$ に対して $x_i$ が互いに独立でかつ同じ確率分布から生起しているとすれば、$D = \{x_1, \cdots, x_N\}$ というデータが手元にあるとして、事後分布

$$ p(\theta | D) = \frac{p(D|\theta)p(\theta)}{p(D)} $$

と書くことが出来ます。 この事後分布を使って、下記のベイズ予測分布を計算し

$$ p(x_{new} | D) = \int_{\theta} p(x_{new} | \theta) p(\theta | D)d\theta $$

真の分布を近似してみたり、あるいは予測に使ってみたりします。これがベイズの基本的な流れになるわけです。

ベイズ予測分布を計算するためには、当然事後分布 $p(\theta | D)$ を予め求めておく必要があり、これを求める（推論する）ことを学習とひとまず呼ぶことにします（実は事後分布を周辺化して消してしまうことを考えれば、予測分布を直接求めるようなアルゴリズムを考えても良いわけで、ベイズ推論は実はもっとフレキシブルです）。

変分ベイズ学習

変分法の心

一旦、ベイズ推論における学習のことに集中しましょう。今達成したいことは、$p(\theta | D)$ という関数を求めるということです。この関数は事後分布と呼ばれてます。この関数を決めることによって$\theta$ の実現値の発生源を決めようとしているわけです（$\theta$を求めようとしているのではなく、 $\theta$ がどのようにバラつくのかを知ろうとしている）。$\theta$ の発生源として使える関数をどのように決定すると良いのでしょうか。

関数を決定する上で大事だと思われるのは、この関数がどんなふうになれば「良い（あるいは悪い）」と言えるかの指標です。そのような指標はすでに与えられているとしましょう。指標を $M(p)$ と書くことにします。仮に良さを表す$M(p)$がどういうものかわかっているなら、これが最大化されるように関数 $p$ を調整してやればいいということです（通常はELBOと呼ばれる指標が使われるが、これは後述する）。

兎にも角にも、$\theta$の発生源となる関数 $p$ を指標 $M(p)$ を使って決めてやろうとしているのです。決して $\theta$ という値をどうするかを決めようとしているのではなく、 $p$ という関数そのものをどのようなものにしてやろうかを決めようとしていることに注意してください。ちなみに、 $M(p)$ は 関数$p$ を引数に取って、スカラー値を返します。このスカラー値が大きければ大きいほど「良い $p$」を入力したということになります。

では関数 $p$ を手当たり次第に試しに入れまくってやって、$M(p)$ を評価していき、一番大きな値を返した $p$ を選べば解決じゃないか！という話です。ただし、当然 $p$ の選び方なんて無数にあるわけなので、この世に存在するありとあらゆる関数（今回の場合は確率分布）を評価して見るなんて無理なのです。

そこで、変分法が登場します。変分法とは、 $M(p)$ に 関数$p$ を試しに入れまくる際に、入れ方を関数の形をほんの少しだけ変えて入れてみるという方法を取ることにします。すなわち、とある $p$ で $M(p)$ を評価したあとに、今度はちょっとだけ $p$ を変えた $ p + \delta p$ で $M(p+\delta p)$ を評価するということです。

もしも上記のことを試してみて、

$$ M(p+\delta p) - M(p) = 0 $$

ということが起こった場合にはどういうことが考えられるでしょう。 $p$ をほんの少しだけ変えても評価値は変わらなかったので、この辺りの $p$ が一先ず解の候補になりそうです。仮に評価指標 $M(p)$ が単峰ならば、上記の $p$ を答えとしてしまうことができます。上記の $M(p+\delta p) - M(p) $ を変分とよんでみることにしましょう。この関数を引数に取ってスカラーの値を返すようなもの（今回で言えば評価指標 $M(p)$ )を「汎関数」と言います。そして、上記の汎関数に入れる関数を少しだけ変えたらどうなるんだろうか？という考えに基づいて変分を考え、関数を決めてやろうという試みを変分法と言います。

なんだか微分と似たような話です！ただし、 $p$ は値ではなく関数なのでもうちょっと話は複雑ですが…。

何が複雑なのかというと、微分であれば、$ f(x) $ の微分を考えるために、$x + \Delta x$ というものを考える必要があります。そのときに $x=5$ に対して $\Delta x = 0.00001$ としてちょっとだけ変えていれてやるか…、みたいなことが出来るわけです。では 関数 $p$ をちょっとだけ変えるために $\delta p$ をどのように選べば良いでしょうか？ 　

　

　 端点を固定しつつ、ちょっとだけ関数の形状をずらす

　 　

　 関数の形状のある一部分のみをずらす 　 上の方が自然な方針に見えますか？それとも下のほうが自然な方針に見えますか？実は私達が求めたいと思っている $p$ の形状の本質を知らない限り、どちらが良いとは言えません（$M(p)$ の評価値を良くするような変更の仕方としか言いようが無いから！）。これが微分に対して変分の方が複雑ですという意味です。（上の方が数学的には扱いやすそうな気はしますよね。物理の１分野である解析力学で学ぶ変分法で最初に学ぶのはおそらくこの形になるでしょう）。

変分ベイズ法の戦略

さて、変分法なるものの存在を知ったところで、変分ベイズ法の話に突入します。まず、私達が求めようとしている関数 $p$ はパラメータ $\theta$ の事後分布です。いま評価指標である $M$ はすでに決まっているとします。上記で見たように、基本的に変分をどのように扱えば良いかの方針というのは簡単には決まりません。$p$ は本来どんな形状であるのかは我々は知りませんし、変分としてどのような変え方をして関数を探していけば良さそうなのかも我々は知りません。

それでも、 $M$ を片っ端から評価していって良さげな $p$ を見つけなければなりません。これを実施するために $p = q(\eta)$ と置いてしまうのです。これは一体何をしたのかというと、 $q(\eta)$ は我々にとって計算しやすく、 $\eta$ を少し変えてやれば $q(\eta)$ の形も少し変わるような良い関数を選んだのです。そうすることで、$p$ の考えうる形状と、そして変分としてどのように変更を加えていくかの方針を $\eta$ に託すことにしたのです。

1. $M(p)$ に基づいて一番良い $p(\theta | D)$ を探したい。
2. $p1(\theta | D), p2(\theta | D),...$ と片っ端から試すのは苦労する。
3. $p$ をほんの少し $\delta p$ だけ変えて、$M(p + \delta p) - M(p)$ を評価していくことに。
4. 関数の探す範囲と変え方 $\delta p$ の方針が分からない。
5. 関数 $p(\theta | D) = q(\eta)$ と置いてしまって、$\eta$ を変えることで変え方と探索する範囲を決定

ということです。そうして得られた $q(\eta)$ をあたかも $\theta$ の事後分布として扱って以後予測分布などに使っていくということになります。この方針で優れている点は、最終的には

$$ M(p(\theta | D)) = M(q(\eta)) $$

と置いてしまったので、変分法で $M(p+\delta p) - M(p) $ を考えることが

$$ \begin{align} & M(p+\delta p) - M(p) \\ \leftrightarrow & M(q(\eta)+\delta q(\eta)) - M(q(\eta)) \\ \leftrightarrow & M(q(\eta + \Delta \eta)) - M(q(\eta)) \\ \leftrightarrow & M(\eta + \Delta \eta) - M(\eta) \end{align} $$

のように書き換えられてしまうことです（$q(\eta)$ は私達が計算しやすく仮定した関数なので、$M(q(\eta))$ という入れ後になった関数は簡単に $M(\eta)$ と入れ子を意識しなくても良い関数に出来てしまうのです）。私達が考えているのは、結局$\eta$ を少し変えることだけなのです。これは従来ニューラルネットワークなどでも使ってきた勾配法が使えます。

そうして $M(\eta)$ を評価して求まった $\eta$ を使って、元々仮定した $q(\eta)$ を $\theta$ の事後分布の形状として採用してやればいいということになります。