分布推定アルゴリズムを学ぶ (2) ～PBIL～

[分布推定アルゴリズムを学ぶ (1) ～概要～](https://qiita.com/_tetsuya_/items/0cc2123f87007e91bd1c)の続き

## Population-Based Incremental Learning(PBIL)とは
EDAsの手法の一つ．
最適化対象は主にビット列であり，変数間の依存性を考慮しない手法．

PBILはベーシックなGenetic Algorithm(GA)とCompetitive Learning(CL)を組み合わせた手法である．そのため，GAの拡張とみなすことができる．基本的にはGAとCLの両方とEDAsのアナロジーを考え，実装に落とし込んでいる．
もともと，EDAsの考え方自体がGAの観点でも説明することができ，EDAsは確率モデル構築型GA(Probabilistic Model-Building GA; PMBGA)と呼ばれることもある．

## 概要

#### Competitive Learning(CL)
CLはいくつかの点を任意個のグループにクラスタリングするために用いられる手法である．各グループに属する点の集合は，任意の類似尺度を用いて決定される．
CLは事前にグループの数を必要としない教師なし学習である．

CLの一般的なネットワークを以下に示す．
入力が各点に対応し，出力が各グループに対応する．このネットワークでは，Inhibitory ConnectionsとExcitatory Connectionsという二種類の接続方式を持つ．
Inhibitory Connectionsは出力ユニット間の接続に用いられ，各世代で活性化するユニットを一つ決定する．
Excitatory Connectionsは，入力から出力への接続部分に用いられる．

![competitive_learning_network.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214163/1b8e3e7c-b9fa-6d1f-101b-5dee445c7742.jpeg)

以下に流れを示す．
まず，初期の重みをランダムに決定する．
つぎに各出力ユニットを以下の計算式によって計算する．

```math
output_i = \sum_j w_{i,j}\times input_j
```
CLでは最も大きな値を持つ出力ユニット(以下，活性化ユニット)のみが入力点に対してアクセスできる．活性化ユニットは入力点の属するグループという意味合いになる．
そして，活性化ユニットに繋がる重みを以下の式によって更新する．

```math
\Delta w_{i,j} = lr\times (input_j - w_{i,j})
```
ここで$lr$は学習率である．

学習はネットワークが安定になるまで繰り返し行われる．
学習が終わると，各出力につながっている重み(重みベクトル)はそのグループを特徴づけるものと解釈できる．
すなわち，重みが大きいという属性は出力ユニットによって表現されたクラスを特徴づける性質になっている．
これがPBILの議論の中心となる概念である．
また，教師ありのCLについても研究がされており，以降では教師なし・教師ありのPBILについてそれぞれ紹介する．


#### 集団から確率モデルへの置き換え
一般的なGAでは候補解の集団を保持し，集団を進化させることで関数空間を探索する．
PBILではその集団を確率ベクトルに置き換えて探索を実現する．正確には，PBILは次世代の集団を作るための候補解を生成するための確率ベクトルを現世代の集団から構築する．

各世代で全ての個体を入れ替え，評価値に比例した選択方法と一般的な交叉を行うGAを仮定した時，$i$番目の変数の値が$j$になる確率は以下で計算される．

```math
P(i,j)=P(\boldsymbol{c}_i=j)=\frac{\sum_{\boldsymbol{v}\in Pop_{G-1}\land v_i=j}EvaluateVector(\boldsymbol{v})}{\sum_{\boldsymbol{v}\in Pop_{G-1}}EvaluateVector(\boldsymbol{v})}
```

ここで$Pop_{G-1}$は$G-1$世代目における集団を意味する．$EvaluateVector$は候補解$\boldsymbol{v}$の評価値を計算する関数である．
この式は評価値による重み付けをした単純な出現割合を求めており，この式によって$(C_{max}-1)\times D$ 確率行列が生成される．ここで$C_{max}$は各変数のとりうる選択肢の数である．$C_{max}-1$になるのは，最後の選択肢の確率は$1-\sum_{k=1}^{C_{max}-1}p_k$で求められるからである．したがって，ビット列($C_{max}=2=|\\{0, 1\\}|$)のときは$D$次元の確率ベクトルとなる．

上記の式を用いて集団を確率行列で表現する際，以下のことに留意する必要がある．

1. 集団としては異なるもの同士でも，確率行列に変換した際には同一のものになる可能性がある．
2. 仮に確率行列を次世代の集団を構築するために用いるとして，生成された候補解が現世代よりも改善しているという保証はない．これは各変数は独立であるという暗黙の了解を置いているにも関わらず，一般の交叉オペレーションではこれを仮定していないからである．

上記の式によって構築した確率行列を次世代の集団を生成するために用いることは不便である(なぜなのかはよく分からない．計算が多少複雑などの理由か？)．代わりにPBILでは以下の表現によって代替する．

集団サイズを$N$としたとき，評価値に基づいた確率的な選択(トーナメント選択など)を用いて$N$個の候補解を集団から復元抽出でリサンプリングする．そして，リサンプリングによって生成された集団を用いて，各変数での共起頻度を用いて確率行列を構築する．

例を以下に示す．
各世代は4個の候補解を持つ集団になっており，各変数(各次元)の共起頻度を用いて確率ベクトルを計算している．
また，Population #1とPopulation #3は異なる集団にも関わらず確率ベクトルの表現上は同じものになっている．これが上記の留意点の一つ目のことである．

![probability_representation.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214163/52cc9867-a29e-2d1d-bfc7-b090578d2d13.jpeg)

PBILはGAとは異なり，集団に対して何かしらの操作(交叉や突然変異)は行わず，確率ベクトルに対してそれらの操作を行う．
PBILは確率ベクトルの各確率を徐々に評価値が良くなる方向に動かしていくという観点においてCLと同様である．
そこで，CLの重みベクトルの更新則をPBILの確率ベクトルの更新則に適用する．すなわち，以下の式で確率ベクトルを更新する．

```math
probability_i=(probability_i\times (1.0-lr))+(lr\times vector_i)
```

ここで$lr$は学習率，$probability_i$は$i$番目の変数が1になる確率，$vector_i$は候補解の$i$番目の値を表している．
この式では，候補解の$i$番目の値が1の場合には確率ベクトルが増加し，0の場合には減少する．

ここまでの議論はGAの交叉に対応する内容であった．
GAの探索において重要なもう一つの操作として突然変異がある．突然変異は集団の多様性を保つために用いられ，特に探索の後半で重要な役割を果たす．なぜなら探索が進むごとに集団内部は同じ候補解が複数存在するようになり，そうなると交叉の効果が薄まるからである．
GAとのアナロジーを考えると，PBILにおいても突然変異の操作を導入することが有用であると考えられる．確率ベクトルに対する突然変異としては以下の二種類の方法がある．

1. 生成された候補解の各変数に対してランダムに値を変更する．
2. 確率ベクトルの各確率に対してランダムに値を増加・減少させる．

PBILでは後者で実装している．

以上をまとめるために，PBILのアルゴリズム全体を以下に示す．
$lr$は学習率，$MUT\\_PROBABILITY$は突然変異確率，$MUT\\_SHIFT$は突然変異時の移動量，$Random(0,1]$は$(0,1]$の範囲の一様乱数，$Random(0~or~1)$は0か1をランダムに出力する乱数生成器である．

![pbil_algorithm.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214163/0d3752f1-973d-7946-54f8-1990a445629b.png)


#### PBILの拡張
以降，ここまでに紹介したPBILをBasic-PBIL(B-PBIL)と呼称する．
B-PBILでは，各世代で複数の個体を生成するにも関わらず評価値が最良の個体一つしか確率ベクトルの更新に利用しない．これはCLの更新則とのアナロジーのためである．
しかし，EDAsの関連研究では最良個体と最悪個体の二つを用いて効率的に確率ベクトルを更新する手法が存在する．
最良個体と最悪個体の各変数を比較することで，どの変数が評価値の良し悪しに影響しているのかという情報を利用できるようになる．そこで，PBILにおいても最良個体と最悪個体を用いた更新則を導入する(これをExtended-PBIL; E-PBILとする)．これによって，B-PBILではどの変数が評価値を良くしているか分からない状況の教師なしCLの更新則だったのに対し，E-PBILではどの変数が評価値を良くしているか分かる(ある種の教師データ)状況の教師ありCLの更新則になる．

E-PBILにおける確率ベクトルの更新則は最良個体に関してはB-PBILと同様である．
しかし，E-PBILではこれに追加して最良個体と最悪個体を比較して値が異なる変数のみ更に確率ベクトルの更新をおこなう．

E-PBILのアルゴリズムを以下に示す．
B-PBILとの主な違いは8～10行目で最良個体と最悪個体の値を比較して異なる場合には$negative\\_lr$の学習率で更に確率ベクトルを更新する部分である．

E-PBILの欠点はユーザーパラメータ($negative\\_lr$)が増えるためにチューニングが面倒な点である．

![extended_pbil_algorithm.PNG](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/214163/a0258562-1006-1519-2ce1-51d0e0300ceb.png)


## 実装
プログラム全体は[github](https://github.com/e5120/EDAs)に置いてある．

### 環境
- windows10
- Python 3.7.3
- anaconda Command line client (version 1.7.2)
- conda 4.7.10
- numpy 1.16.2

[github](https://github.com/e5120/EDAs)の方に[Singularity](https://sylabs.io/guides/3.2/user-guide/index.html)の環境ファイルを置いてあるので，それをビルドしても動く(動作確認済み 2019/12/19時点)．

### コード
今後PBIL以外の手法も実装していくので，EDAs用の基底クラスを以下に示す．
コンストラクタの引数`categories`は目的関数の次元サイズのリストであり，各要素は対応する次元の選択肢の数を表現している．
例えば，5次元ビット列の場合には，`categories=[2,2,2,2,2]`となる．
また，`lam`は集団サイズを指し，`theta_init`は初期の確率分布を明示的に指定するための変数である(基本はNone，このとき一様分布)．

`sampling`メソッドでは，確率モデルから個体を1つサンプルする．個体はone-hotベクトルの形式である．

```python:eda_base.py
import numpy as np


class EDABase(object):
    def __init__(self, categories, lam=2, theta_init=None):
        self.N = np.sum(categories - 1)

        self.d = len(categories)
        self.C = categories
        self.Cmax = np.max(categories)
        self.theta = np.zeros((self.d, self.Cmax))
        for i in range(self.d):
            self.theta[i, :self.C[i]] = 1.0 / self.C[i]
        self.valid_params = int(np.sum(self.C - 1))
        self.valid_d = len(self.C[self.C > 1])

        if theta_init is not None:
            self.theta = theta_init

        self.lam = lam
        # for record
        self.best_eval = np.inf
        self.best_indiv = None
        self.eval_count = 0

    def sampling(self):
        rand = np.random.rand(self.d, 1)
        cum_theta = self.theta.cumsum(axis=1)

        c = (cum_theta - self.theta <= rand) & (rand < cum_theta)
        return c

    def is_convergence(self, eps=1e-8):
        return np.abs(self.theta.max(axis=1).mean() - 1.0) < eps

    def update(self, c_one, fxc, range_restriction=False):
        raise NotImplementedError
```

PBILの実装を以下に示す．
基本的には擬似コードの通りだが，`update`メソッドの最後の8行のみが異なる．
これは，確率分布が0，または1に完全に収束しないように目的関数の次元に応じてクリッピングする処理である．`range_restriction=True`の場合のみクリッピングされる．
論文の擬似コードでは，学習率の値とサンプリングされた個体の偏りによっては，特定の変数の確率分布が誤った値に早期収束する可能性がある(と思っている)．
そのため，独自にクリッピング処理を追加している．これに関する考察は後述の実験で行う．

また，`sampling`メソッドで個体をone-hot形式で返しているため，`c_one`は`(集団サイズ,次元数,one-hot)`の形のndarrayを想定している．

```python:pbil.py
import numpy as np

from optimizer.eda_base import EDABase


class PBIL(EDABase):
    def __init__(self, categories, lr, lam=64, negative_lr=None,
                 mut_prob=0.02, mut_shift=0.05, theta_init=None):
        super(PBIL, self).__init__(categories, lam=lam, theta_init=theta_init)
        # only binary strings
        assert self.Cmax == 2
        assert 0.0 < lr < 1.0
        assert negative_lr is None or 0.0 < negative_lr < 1.0
        assert 0.0 <= mut_prob <= 1.0
        
        self.lr = lr
        self.negative_lr = negative_lr
        self.mut_prob = mut_prob
        self.mut_shift = mut_shift

    def update(self, c_one, fxc, range_restriction=False):
        self.eval_count += c_one.shape[0]
        # sort by fitness
        idx = np.argsort(fxc)
        best_idx = idx[0]
        # store best individual and evaluation value
        if self.best_eval > fxc[best_idx]:
            self.best_eval = fxc[best_idx]
            self.best_indiv = c_one[best_idx]
        # update probability vector
        best_indiv = c_one[best_idx]
        self.theta[:, -1] = (1.0 - self.lr) * self.theta[:, -1] + self.lr * best_indiv[:, -1]
        # update probability vector by negative sample
        if self.negative_lr is not None:
            worst_idx = idx[-1]
            worst_indiv = c_one[worst_idx]
            diff_dix = best_indiv[:, -1] != worst_indiv[:, -1]
            self.theta[diff_dix, -1] = (1.0 - self.negative_lr) * self.theta[diff_dix, -1] + self.negative_lr * best_indiv[diff_dix, -1]
        # mutation probability vector
        dim = self.theta.shape[0]
        mut_idx = np.random.rand(dim) < self.mut_prob
        mut_num = np.sum(mut_idx)
        self.theta[mut_idx, -1] = (1.0 - self.mut_shift) * self.theta[mut_idx, -1] \
                                  + np.random.randint(0, 2, mut_num) * self.mut_shift
        self.theta[:, 0] = 1 - self.theta[:, -1]
        # if range_restriction is True, clipping
        for i in range(self.d):
            ci = self.C[i]
            theta_min = 1.0 / (self.valid_d * (ci - 1)) if range_restriction and ci > 1 else 0.0
            self.theta[i, :ci] = np.maximum(self.theta[i, :ci], theta_min)
            theta_sum = self.theta[i, :ci].sum()
            tmp = theta_sum - theta_min * ci
            self.theta[i, :ci] -= (theta_sum - 1.0) * (self.theta[i, :ci] - theta_min) / tmp
            self.theta[i, :ci] /= self.theta[i, :ci].sum()
```

## 目的関数
ここではEDAsの性能評価においてよく用いられる三つのベンチマーク関数について紹介する．

#### One-Max問題
以下で定義される目的関数を指す．

```math
f(\boldsymbol{c})=\sum_{i=1}^{D}c_i \\
```

ここで，$\boldsymbol{c}=(c_1,c_2,...,c_D)\in \\{0,1\\}^D$の$D$次元ビット列である．
すなわち，ビット列中の1の総数が評価値になる．
大域的最適解は$\boldsymbol{c}^*=(1,...,1)$である．

#### Four Peaks問題
以下で定義される目的関数を指す．

```math
\begin{align}
f(\boldsymbol{c})&=\max(o(\boldsymbol{c}),z(\boldsymbol{c}))+REWARD \\
o(\boldsymbol{c})&=先頭から連続で1が並んでる数 \\
z(\boldsymbol{c})&=後尾から連続で0が並んでる数 \\
REWARD&=\left\{
\begin{array}{ll}
D & o(\boldsymbol{c}) > T \ \ and \ \ z(\boldsymbol{c}) > T \\
0 & otherwise
\end{array}
\right.
\end{align}
```

ここで，$\boldsymbol{c}=(c_1,c_2,...,c_D)\in \\{0,1\\}^D$の$D$次元ビット列，$T$はしきい値である．

以下に例を示す．
次元 $D=10$とし，しきい値 $T=2$，個体 $\boldsymbol{c}=(1,1,0,0,1,0,1,0,0,0)$の場合．
$o(\boldsymbol{c})=2,\ z(\boldsymbol{c})=3$であり，$o(\boldsymbol{c})>T$を満たさないので，$REWARD=0$となる．
したがって評価値 $f(\boldsymbol{c})=\max(2,3)+0=3$となる．

この例の場合，大域的最適解は以下の二つである．

```math
\boldsymbol{c}_1^*= (1,1,1,0,0,0,0,0,0,0) \\
\boldsymbol{c}_2^*= (1,1,1,1,1,1,1,0,0,0)
```

また，局所最適解は以下の二つである．

```math
\boldsymbol{c}_1= (0,0,0,0,0,0,0,0,0,0) \\
\boldsymbol{c}_2= (1,1,1,1,1,1,1,1,1,1)
```

#### 3-Deceptive問題
以下で定義される目的関数を指す．

```math
\begin{align}
f(\boldsymbol{c})&=\sum_{i=0}^{\lceil D/3 \rceil -1}g(d,c_{3i+1},c_{3i+2},c_{3i+3}) \\
g(d,c_1,c_2,c_3)&=\left\{
\begin{array}{ll}
1-d & \sum_{j=1}^3c_j=0 \\
1-2d & \sum_{j=1}^3c_j=1 \\
0 & \sum_{j=1}^3c_j=2 \\
1 & \sum_{j=1}^3c_j=3 \\
\end{array}
\right.
\end{align}
```

ここで，$\boldsymbol{c}=(c_1,c_2,...,c_D)\in \\{0,1\\}^D$の$D$次元ビット列，$d$はユーザーパラメータで一般に$d=0.1$が用いられる．

以下に例を示す．
次元 $D=9$とし，$d=0.1$，個体 $\boldsymbol{c}=(1,1,0,0,1,0,0,0,0)$の場合．
まず，以下のように三つの部分列に分割される．

```math
\begin{align}
(c_1,c_2,c_3)&=(1,1,0) \\
(c_4,c_5,c_6)&=(0,1,0) \\
(c_7,c_8,c_9)&=(0,0,0) \\
\end{align}
```

各部分列を関数$g$に入力すると，

```math
\begin{align}
g(d,c_1,c_2,c_3)&=0 \\
g(d,c_4,c_5,c_6)&=0.8 \\
g(d,c_7,c_8,c_9)&=0.9 \\
\end{align}
```

となり，評価値 $f(\boldsymbol{c})=0+0.8+0.9=1.7$と求められる．

この例の場合，大域的最適解は$\boldsymbol{c}^*=(1,1,1,1,1,1,1,1,1)$，局所最適解は1つ以上の部分列が$(0,0,0)$となる場合である．


## 実験(1)：B-PBILの性能評価
### 設定
目的関数はOne-Max問題，Four Peaks問題，3-Deceptive問題を用いる．
One-Max問題は$D=100$，Four Peaks問題は$D=\\{50,100\\}$，3-Deceptive問題は$D=21$とする．
また，Four Peaks問題のしきい値 $T=0.1D$，3-Deceptive問題のユーザーパラメータ $d=0.1$とする．

その他のユーザーパラメータは以下の通りである．

|パラメータ名 |値 |
|---|---|
|集団サイズ：$N$ |{8, 32, 64, 128} |
|学習率：$lr$ |0.1 |
|クリッピング処理 | False |
|突然変異確率：$MUT\\_PROBABILITY$ |0.02 |
|突然変異時の移動量：$MUT\\_SHIFT$ |0.05 |

以下のいずれかを満たした場合，最適化を終了する．
1. 大域的最適解をサンプルできたとき(成功)
2. 確率分布が収束したとき(失敗)
3. 上限世代数($2\times 10^3$)に達したとき(失敗)

以上の設定で，乱数のシードを変更して10試行行う．

### 観察指標
以下の指標を記録する．

1. 成功試行数(大域的最適解を見つけられた回数)
2. 平均評価回数($\pm$標準偏差)
3. 最適化終了時の最良評価値の平均($\pm$標準偏差)

ここで平均評価回数は，成功試行時のみのものとする．
評価回数は以下で計算する．

```math
評価回数 = 集団サイズ \times 世代数
```

### 結果
One-Max問題($D=100$)をB-PBILで最適化したときの結果を以下の表に示す．

|集団サイズ |成功試行数 | 平均評価回数($\pm$標準偏差)|平均最良評価値($\pm$標準偏差)|
|---|---|---|---|
|8 |10/10 |1566.4$\pm$882.40|100.0$\pm$0.0|
|32 |10/10 |2268.8$\pm$138.34|100.0$\pm$0.0|
|64 |10/10 |3840.0$\pm$458.84|100.0$\pm$0.0|
|128 |10/10 |6579.2$\pm$632.79|100.0$\pm$0.0|

Four Peaks問題($D=50$)をB-PBILで最適化したときの結果を以下の表に示す．

|集団サイズ |成功試行数 | 平均評価回数($\pm$標準偏差)|平均最良評価値($\pm$標準偏差)|
|---|---|---|---|
|8 |10/10 |7176.0$\pm$2633.44|94.0$\pm$0.0|
|32 |10/10 |8912.0$\pm$6083.41|94.0$\pm$0.0|
|64 |10/10 |14924.8$\pm$13962.27|94.0$\pm$0.0|
|128 |10/10 |9984.0$\pm$3601.78|94.0$\pm$0.0|

Four Peaks問題($D=100$)をB-PBILで最適化したときの結果を以下の表に示す．

|集団サイズ |成功試行数 | 平均評価回数($\pm$標準偏差)|平均最良評価値($\pm$標準偏差)|
|---|---|---|---|
|8 |1/10 |13776.0$\pm$0.0|100.9$\pm$53.30|
|32 |4/10 |38192.0$\pm$11244.68|129.8$\pm$49.04|
|64 |2/10 |89568.0$\pm$1120.00|122.2$\pm$35.28|
|128 |5/10 |111539.2$\pm$44034.13|144.0$\pm$45.00|

3-Deceptive問題($D=21$)をB-PBILで最適化したときの結果を以下の表に示す．

|集団サイズ |成功試行数 | 平均評価回数($\pm$標準偏差)|平均最良評価値($\pm$標準偏差)|
|---|---|---|---|
|8 |0/10 |-|6.49$\pm$0.12|
|32 |0/10 |-|6.56$\pm$0.09|
|64 |0/10 |-|6.73$\pm$0.09|
|128 |0/10 |-|6.73$\pm$0.08|


### 考察
One-Max問題の場合，集団サイズが小さい方が評価回数が少なくて済むという結果になった．
これは，大量にサンプルしても更新に用いられる個体は一つのみであり，それならば少ない個体で多くの世代を回した方が確率分布の収束に合わせながら更新できるので有利なためと考えられる．
しかし$N=8$の場合，更新に用いられる最良個体が不安定になってしまい，標準偏差が大きくなってしまっている．
そのため，このユーザーパラメータ群の場合は評価回数と安定性の観点で$N=32$程度が適切と考えられる．

Four Peaks問題では，$D=50$の場合にはどの集団サイズでも全試行成功しているが，$D=100$の場合は最大でも5試行しか成功していない．
これは，評価値の$REWARD=D$となるための条件が難しくなり，局所最適解に収束している試行があるためである．
One-Max問題の場合はとりあえず1になる変数を増やせば良かったが，Four Peaks問題の場合は個体の先頭は1にしつつ，同時に後方は0にしなければならない．そのような学習のためには$REWARD$の条件をクリアする個体が集団内に一つでも存在しつづけなければならない．
そのため，一世代で多くの個体をサンプルできる集団サイズが大きい方が成功試行が多い傾向にあると考えられる．

3-Deceptive問題の場合，どの集団サイズでも全試行失敗している．
これもFour Peaks問題($D=100$)と同じような理由であり，3-Deceptive問題の場合騙し構造になっている部分問題を解かなければならず，各変数を独立に学習しているPBILでは解くことのできない問題設定であると考えられる．
ただし，集団サイズが大きいほうが部分列$(1,1,1)$を初期段階でサンプルできる確率が上がるため，平均最良評価値は集団サイズの増加にしたがって改善傾向にある．

## 実験(2)：B-PBILとE-PBILの性能比較

### 設定
目的関数は$D=100$のOne-Max問題を用いる．

その他のユーザーパラメータは以下の通りである．

|パラメータ名 |値 |
|---|---|
|集団サイズ：$N$ |32 |
|学習率：$lr$ |0.1 |
|負の学習率：$negative\\_lr$(E-PBILのみ) |{0.01, 0.075, 0.1, 0.2} |
|クリッピング処理 | False |
|突然変異確率：$MUT\\_PROBABILITY$ |0.02 |
|突然変異時の移動量：$MUT\\_SHIFT$ |0.05 |

終了条件は実験(1)と同様

### 観察指標
実験(1)から平均最良評価値を除いたもの．

### 結果
One-Max問題($D=100$)をB-PBILとE-PBILで最適化したときの結果を以下の表に示す．

|手法 |$negative\\_lr$ |成功試行数 | 平均評価回数($\pm$標準偏差)|
|---|---|---|---|---|
|B-PBIL |- |10/10 |2268.8$\pm$138.34|
|E-PBIL |0.01 |10/10 |2710.4$\pm$1638.52|
|E-PBIL |0.075 |10/10 |2438.40$\pm$1671.78|
|E-PBIL |0.1 |10/10 |2428.8$\pm$1707.46|
|E-PBIL |0.2 |10/10 |1785.6$\pm$817.82|

### 考察
結果より，基本的にB-PBILが最もいい結果になっている．E-PBILの$negative\\_lr=0.2$が平均最良評価値の観点では最もいい結果になっているが，標準偏差まで考慮するとB-PBILとそこまで変わらない．
E-PBILで評価回数が改善しなかった理由は，学習初期に間違った更新の場合にも大きく更新されてしまうためと考えられる．
つまり，ある変数において最良個体では0で最悪個体では1の時に確率ベクトルが0の方に二回移動するために，初期の学習が遅れてしまっているためである．
学習初期はほとんど一様分布からのサンプリングのため，今回の設定($N=32$)の場合には最良個体としてあまり評価値の良くない個体が選ばれてしまう場合がある．
集団サイズを大きくすれば，各世代での最良個体がある程度良いものになると考えられるのでE-PBILの方が評価回数が改善する可能性がある．
今回は実験(1)より，One-Max問題の場合$N=32$程度が良さそうであると事前に分かっていたのでその設定で実験を行った．

## 実験(3)：突然変異の影響
### 設定
目的関数は$D=100$のOne-Max問題を用いる．

その他のユーザーパラメータは以下の通りである．

|パラメータ名 |値 |
|---|---|
|集団サイズ：$N$ |32 |
|学習率：$lr$ |0.1 |
|クリッピング処理 | False |
|突然変異確率：$MUT\\_PROBABILITY$ |{0.0, 0.02, 0.1} |
|突然変異時の移動量：$MUT\\_SHIFT$ |{0.01, 0.05, 0.1} |

終了条件は実験(1)と同様

### 観察指標
実験(1)から平均最良評価値を除いたもの．

### 結果
One-Max問題($D=100$)をB-PBILで最適化したときの結果を以下の表に示す．
$MUT\\_P$は$MUT\\_PROBABILITY$，$MUT\\_S$は$MUT\\_SHIFT$を意味している．

|$MUT\\_P$ |$MUT\\_S$|成功試行数 | 平均評価回数($\pm$標準偏差)|
|---|---|---|---|
|0.0 | -| 7/10 |2057.14$\pm$94.80|99.60$\pm$0.66|
|0.02 | 0.01| 10/10 |3369.6$\pm$2663.28|100.0$\pm$0.0|
|0.02 | 0.05| 10/10 |2268.8$\pm$138.34|100.0$\pm$0.0|
|0.02 | 0.1| 10/10 |2371.2$\pm$223.75|100.0$\pm$0.0|
|0.1 | 0.01| 10/10 |2272.0$\pm$396.34|100.0$\pm$0.0|
|0.1 | 0.05| 10/10 |2720.0$\pm$285.50|100.0$\pm$0.0|
|0.1 | 0.1| 10/10 |5302.4$\pm$2243.52|100.0$\pm$0.0|


### 考察
突然変異を用いない，つまり$MUT\\_PROBABILITY=0.0$の場合3試行失敗しており，突然変異の有効性が確認された．
突然変異を用いる場合は，確率と移動量を同時に極端に大きくすると評価回数の観点で大きく劣化した．
一般のGAと同様に突然変異の影響が大きすぎると，上手く学習できている部分が破壊されてしまい悪影響のほうが大きくなってしまうためである．

## 実験(4)：学習率と確率ベクトルのクリッピングの影響
### 設定
目的関数は$D=100$のOne-Max問題を用いる．

その他のユーザーパラメータは以下の通りである．

|パラメータ名 |値 |
|---|---|
|集団サイズ：$N$ |32 |
|学習率：$lr$ |{0.01, 0.1, 0.3, 0.5} |
|クリッピング処理 | {False, True} |
|突然変異確率：$MUT\\_PROBABILITY$ |0.02 |
|突然変異時の移動量：$MUT\\_SHIFT$ |0.05 |

終了条件は実験(1)と同様

### 観察指標
実験(1)から平均最良評価値を除いたもの．

### 結果
One-Max問題($D=100$)をB-PBILで最適化したときの結果を以下の表に示す．

|$lr$ |クリッピング|成功試行数 | 平均評価回数($\pm$標準偏差)|
|---|---|---|---|
|0.01 | False| 10/10 |27993.6$\pm$2795.41|$\pm$|
|0.01 | True| 10/10 |28694.4$\pm$3982.43|
|0.1 | False| 10/10 |2268.8$\pm$138.34|
|0.1 | True| 10/10 |2355.2$\pm$449.19|
|0.3 | False| 10/10 |11363.2$\pm$5515.27|
|0.3 | True| 10/10 |1843.2$\pm$297.51|
|0.5 | False| 10/10 |12419.2$\pm$3473.19|
|0.5 | True| 10/10 |1910.4$\pm$412.80|


### 考察
クリッピング処理を用いない場合，$lr=0.1$が最も良い結果となっている．
一般に集団サイズと学習率は反比例の関係にあることが知られており，$N=32$の場合には$lr=0.1$程度がちょうどいい値であったのだろう．
$lr=0.01$の場合には一回の更新量が少なすぎて収束に時間がかかってしまい，$lr=\\{0.3,0.5\\}$の場合には一回の更新量が大きすぎて誤った方向に更新した際の悪影響が大きいために評価回数が増えていると考えられる．
一方で，クリッピング処理を適用した際には$lr=0.01$を除いて平均評価回数が2000程度で安定した．
これは，誤った方向に更新されて間違った値に確率が収束しても，クリッピングによって正しい値がサンプルされる確率もある程度残るために誤った更新の悪影響が緩和されているためと考えられる．
確率ベクトルのクリッピング処理は多様性を強制的に保存しているとみなすことができ，初期に誤った方向に学習した変数も終盤に正しい方向に学習し直しやすかったと考えられる．

## まとめ
PBILの概要について説明し，自分の気になった部分に関して実験を行った．
実験より以下の知見を得た．

1. PBILは(One-Max問題のような)単純な問題設定の場合には非常に有効な手段．
2. 今回の設定ではB-PBILとE-PBILの間に大きな差はなかった．
3. EDAsにおいても突然変異は多様性を保つための有力な手段である．
4. 誤った方向への確率分布の初期収束を防ぐためにクリッピング処理を独自に導入し，一定の効果を確認した．

## 参考文献
1. [Population-Based Incremental Learning: A Method for Integrating Genetic Search Based Function Optimization and Competitive Learning](https://www.ri.cmu.edu/pub_files/pub1/baluja_shumeet_1994_2/baluja_shumeet_1994_2.pdf)
2. [Removing the Genetics from the Standard Genetic Algorithm](https://www.sciencedirect.com/science/article/pii/B9781558603776500141)
3. [From recombination of genes to the estimation of distributions I Binary parameters](http://www.muehlenbein.org/estbin96.pdf)
4. 遺伝的アルゴリズム -その理論と先端的手法-