この記事は前半が解きたい問題の説明、後半が自分なりに考えた解法の説明です。

問題

0から1の範囲に散らばっているデータをヒストグラムにすると次のような感じになったとします。これは0から1までまんべんなく散らばっています。まんべんなく散らばっているのをこの記事では正常ということにします。

次のヒストグラムでは0.5付近に異常があります。

次のヒストグラムでは0.7付近にピークがありますが、もともとの数が少ないので、偶然ともいえそうです。これは正常ということにします。

次のヒストグラムでは、0.6付近と0.9付近の2箇所に異常があります。

こんな感じで、数値の集合のデータから、異常かどうかの判定に加え、異常な個所はどこなのかを検知したいという場合に、どうしたらいいかを模索しました。

もうちょっと詳しく

次のグラフのオレンジの線で表すの分布を推測して、オレンジの線のピークの箇所を検知します。

与えられるデータは、母集団からたまたま抽出された標本だとみなし、オレンジの線は母集団の分布、青いヒストグラムは標本の分布を表します。標本である数値の集合からもとの母集団の分布を求めることになります。

母集団の分布は一様分布と0個以上の正規分布の重ね合わせとします。

上のグラフの例では、母集団の分布は次の3つの重ね合わせです。

一様分布 79%
平均0.56、標準偏差0.009の正規分布 10%
平均0.89、標準偏差0.010の正規分布 11%

79%などと書いてあるのは分布の重みです。数値がもし100個あれば、一様分布に従う数値は79個、1つ目の正規分布に従う数値は10個、2つ目の正規分布に従う数値は残り11個であることを表します。

次のヒストグラムはピークがない例です。一様分布のみなのでオレンジの線は直線です。とがっているところが異常なピークなのかランダムな振れ幅なのかの区別が難しいのですが、ヒストグラム全体を人間が見て異常にデータが偏っていると即座に感じるもの以外は、ピークがなく異常なしとします。（それでも区別が難しいものはありますが）

問題設定

0.0から1.0までの数値の集合である標本から、母集団の分布を求める。母集団の分布は0.0から1.0までの一様分布と0個以上の正規分布の重ね合わせとする。

ヒストグラムを人間が見て判別できる程度の明瞭な分布のみを検知できればよいものとする。以下の検知漏れは許容か、むしろ検知しないことを推奨する。

ヒストグラムを人間が見ても、すべての正規分布のピークを識別できないほどに正規分布の数が多すぎるケースでは、一部のピークのみを検知するか、まったく検知せずに一様分布のみと判定してもやむなし
ヒストグラムを人間が見ても、ピークが正規分布によるものなのか、一様分布の中のランダムな振れなのかが判別できないほどに、ピークの高さが低いケースでは、正規分布を検知しなくてもやむなし
ヒストグラムを人間が見ても、明瞭なピークを認識できないほど広がりを持った（分散の大きい）正規分布は、検知しなくてもやむなし

入力は標本のデータで、数値の集合。

出力は母集団の分布を表すのに必要なパラメータで、次の通り。

一様分布の重み
正規分布の数
正規分布ごとに
- 平均
- 標準偏差
- 重み

アルゴリズム

ここからは考えたアルゴリズムの説明です。

解（母集団の分布）の候補に対して標本と一致する度合いをスコアとして計算できるロジックを決め、そのスコアが最大となるような解候補を探します。

不自然な解を最適としないようスコアリングのロジックが鍵になります。いろいろと試行錯誤の末、以下の解法のスコアリングが私としては一番納得できるものになりました。

最適解の探索アルゴリズムとしては、焼きなまし法とビームサーチを組み合わせたようなイメージですが、この部分はこれがいい方法なのかはよくわかりません。工夫の余地はたくさんあると思います。

※一様分布がなく正規分布のみであれば、k-meansやEMアルゴリズムによる混合ガウスモデルのクラスタリングでできます。この問題は一様分布も混ざっている場合にどうしたらいいか、ということになります。

全体ループ

分布を表すパラメータをランダムに決めた母集団分布候補を複数作成する
母集団分布候補ごとに、パラメータをランダムに少しだけ変えた複数の新しい母集団分布候補を作成する
新しい母集団分布候補ごとに、そのスコアを計算する
新しい母集団分布候補群から一定数を選び、それを母集団分布候補群として、2に戻る。一定回数繰り返したら終了する

候補群から一定数を選択する方法

母集団分布候補をスコア順に並べる
スコアのもっともよい候補をまずは選択する
ランダムに残りを選択する。その際はスコアのよい候補を優先する

スコア順に並べて上位から機械的に一定数を選択するのではなく、ランダムに選択するようにします。ランダム性は上位ほど選ばれやすいようにします。つまり上位は選ばれる確率は高いがたまに選ばれないこともあり、下位は選ばれないことが多いがたまに選ばれるようにします。

全体ループの回り始めは下位も多く選ばれやすいのですが、ループの回数を重ねるにつれ上位に絞られるようになり、後半はほとんど最上位のみしか選ばれなくなります。

候補群はグループに分けて、グループ内でソートして選択し、グループ間でスコアに差があってもグループ単位で絶滅しないようにします。最初はスコアのよいグループもすぐに頭打ちになって、あとからスコアが伸びるグループが追い越すことがあるためです。

ループを繰り返すうちにスコアのいいものが少しずつ見つかっていく様子の例が下のグラフです。縦軸はスコア、横軸は進捗です。横軸の数値は解候補の番号で、ループ回数より大きな値です。グループで色分けしています。

この例ではグループ間でスコアの伸びの速度が違います。緑のグループが最適解ですが、グルーピングしないと、紫が他を最初に駆逐してしまう可能性があります。そうでなくても次は紺が緑を駆逐してしまいます。

グループは正規分布の数0個から4個までの5グループにわけました。全体ループの中でパラメータをランダムに変えるときには、正規分布の数以外のパラメータを変えます。

このグラフでいうと紫は正規分布なしのグループです。下記スコア計算のところで説明の通り正規分布なしではスコアが常に0になるため、横一直線に並んでいます。紫が最適解なら異常なしという判定になります。

スコアの計算方法

分布の確率密度関数に各数値を代入して、各数値の確率密度を計算する
各確率密度を $p$ として $z=\frac{p-1}{p+1}$ を計算する
$z$ の平均を計算する（算術平均）
その値を $\bar{z}$ 、データ数を $n$ 、自由度を $d$ として $\bar{z}-\frac{kd}{n}$ を計算し、スコアとする。 $k$ は定数でこのアルゴリズムでのハイパーパラメータとなる

スコア計算手順1：確率密度

分布から確率密度関数がわかるので、数値ごとに確率密度を計算します。

次の分布の場合、

一様分布、重み $w_0$
平均 $\mu_1$ 、標準偏差 $\sigma_1$ の正規分布、重み $w_1$
平均 $\mu_2$ 、標準偏差 $\sigma_2$ の正規分布、重み $w_2$

確率密度関数は次のとおりです。

p = w_0 + \frac{1}{\sqrt{2\pi}\sigma_1}e^{-\frac{(x-\mu_1)^2}{2\sigma_1^2}} + \frac{1}{\sqrt{2\pi}\sigma_2}e^{-\frac{(x-\mu_2)^2}{2\sigma_2^2}}

スコア計算手順2：値域変換

手順1で計算した各確率密度を $p$ として $z=\frac{p-1}{p+1}$ を計算します。

確率密度は0より大きな値です。正規分布のピークの部分は1を超え、標準偏差を小さくすればいくらでもピークでの確率密度を大きくすることができます。確率密度の平均をスコアとすると、最適なスコアとなる分布は1個のデータだけにマッチした標準偏差を極端に小さくした正規分布を含むものになってしまいます。これは欲しい結果ではありませんので、これを防ぐために $z=\frac{p-1}{p+1}$ の計算式で補正しています。

値域 $p>0$ を値域 $z>0$ かつ $z<1$ に変換する式として、最初は $\frac{p}{p+1}$ を考えました。これは以下の性質があります。

$p=0$ で $\frac{p}{p+1}=0$
$p$ が0付近では $\frac{p}{p+1}$ は $p$ で近似できる直線に載る
$p$ が大きくなるにつれ $\frac{p}{p+1}$ の傾きがなだらかになり、 $p \rightarrow +\infty$ で $\frac{p}{p+1} \rightarrow 1$

この $\frac{p}{p+1}$ でも十分なのですが、正規分布がなく一様分布のみの場合データの内容に依らず常に $p=1$ となり、 $\frac{p}{p+1}=\frac{1}{2}$ という定数になります。せっかくなのでスコアとしてはこれを0にするために、 $\frac{1}{2}$ 引いて、上限を1のままにするために2倍しました。

z = 2(\frac{p}{p+1} - \frac{1}{2}) = \frac{2p}{p+1} - 1 = \frac{p-1}{p+1}

スコア計算手順3：平均

手順2で計算した $z$ の算術平均を計算します。

尤度との類推からすべてを掛け算する相乗平均または対数の算術平均を最初は試したのですが、あまりいい結果にはなりませんでした。ヒストグラムを人間が見て直感的に判別した分布と違う結果を出すことが、単純な算術平均よりも多かったです。

スコア計算手順4：ペナルティ項

手順3で計算した平均値を $\bar{z}$ 、データ数を $n$ 、自由度を $d$ として $\bar{z}-\frac{kd}{n}$ を計算し、それをスコアとします。 $k$ は定数でこのアルゴリズムにおけるハイパーパラメータとなります。

自由度は分布に含まれるパラメータの数で、一様分布のみであれば自由度0、正規分布が1個なら3、正規分布が2個なら6です。自由度が多いほど複雑な分布を表現でき、スコアも上げることができます。正規分布最大4個の範囲で最適解を求めると、このペナルティ項がない場合、たいてい正規分布4個が最適解となってしまいます。

ヒストグラムを人間が見て直感的に判別した分布というのは、自由度のできるだけ少ないシンプルな分布になりますので、それに寄せるためのペナルティ項が $-\frac{kd}{n}$ です。

自由度が1つ増えるごとに一定の個数のデータの $z$ 値を一定量増加させることができると考えました。「一定の個数×一定量」は、全体のデータ数や分布の内容には依存しないのではないかとも考えました。この「一定の個数×一定数」を $k$ としたとき、自由度 $d$ では $kd$ となり、 $\bar{z}$ はデータ数 $n$ で除算した値なので、 $\bar{z}$ に対するペナルティとしては $\frac{kd}{n}$ となります。 $k$ の値は試行錯誤した結果だいたい $1.0$ から $3.0$ の間あたりがよさそうです。ただ、ひょっとしたらデータ数をずっとずっと大きくした場合には適正値が変わってくるかもしれません。

問題点

k-meansや混合ガウスモデルのクラスタリングに比べて計算コストが桁違いに大きいです。

もっといい方法がないかな。。。

と1か月ぐらい前は思っていました。この記事は書いてから未公開のまま1か月経過したものですが、その間にこの問題から離れて、いま改めて見ると、ふつうに勾配降下法で最適解を求めればいいのではと思ってます。。。せっかくやったので記録の意味でQiitaに残しておきます。