混合ガウスモデルの理論と実装

混合ガウスモデルの理論と実装 scikit-learnを0から作ろうProject 第7回勉強会 (2024-06-28) 小宮和真 武蔵野大学データサイエンス学科 4年 GitHub: misya11p Twitter: @ds33miya

• https://github.com/misya11p
• https://twitter.com/ds33miya

前置き 今回は混合ガウスモデル（GMM: Gaussian Mixture Model）を実装する。このモデルはクラスタリングなどに応用できる。 このモデルを理解するにあたってベイズ統計の知識が必要になる。一応、必要な部分はざっくりと本資料にもまとめているが、おそらく不 十分であるため、分からない点は調べるなり質問するなりしていただければ。 本資料は、まずはじめに混合ガウスモデルの実装に必要な知識を説明し、実際に実装を行う。その後、発展的な内容として、より一般的な 学習アルゴリズムであるEMアルゴリズムを説明する。この発展的な内容は、VAEや拡散モデルといった生成モデルの基礎でもあるため、 その辺りに興味がある人はぜひ学んで欲しい。 2

目次 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 3

ガウス分布 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 4

ガウス分布 確率分布 ある確率変数が取りうる値とその確率を表したもの。 *確率変数: 値が確率的に決まる変数 例) サイコロを1回振った時に出る目を確率変数xとするとき、xの確率分布p(x)は次のようになる。 x 1 2 3 4 5 6 p(x) 1/6 1/6 1/6 1/6 1/6 1/6 確率分布を表す関数p(x)は確率密度関数と呼ぶ。この例だとp(x) = 1/6。 5

ガウス分布 確率分布 確率分布p(x)は観測値xの確率密度を返す関数である。返すものは確率ではないことに注意。ただ確率変数が離散的な場合、確率密度と確 率が一致するため、実質的に確率を返す関数と見られる。 確率変数が連続的である場合、確率変数が取りうる全ての値に確率を定義することができない。そこで、観測値xと範囲[a, b]を与えたと き、xが[a, b]から得られる確率が[a, b]での積分で表せるように確率密度関数を定義する。 b ∫ p(x)dx ​ ​ a 6

ガウス分布 ガウス分布 確率分布の一種。 2 p(x; μ, σ ) = (x − μ)2 ) exp (− 2σ 2 2πσ 2 1 ​ ​ ​ ​ 正規分布と呼ぶことの方が多いかも。ただ混合と付く場合はガウス分布と呼ぶことの方が多い気がするので、本資料でもそう呼ぶ。 また、ガウス分布であることを明示するためにpではなくN やN を使うことも多い（Normal DistricutionのN）。本資料ではN を使う（N はデータ数として使う）。 N (x; μ, σ 2 ) ​ 7

ガウス分布 多変量ガウス分布 確率変数を多次元のベクトルx ∈ Rd としたもの。 p(x; μ, Σ) = 1 (2π)d ∣Σ∣ ​ ​ exp ((x − μ)T Σ−1 (x − μ)) ​ μ ∈ Rd : 平均ベクトル Σ ∈ Rd×d : 共分散行列 8

ガウス分布のパラメータ推定 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 9

ガウス分布のパラメータ推定 確率分布のパラメータ 確率分布の詳細な形を決める変数をパラメータと呼ぶ。 例えばガウス分布は平均μと分散σ 2 の二つのパラメータを持ち、そ れらを変えることで形状が変わる。 パラメータは確率変数の右にセミコロンを挟んで示すことが多い。またパラメータはθ で表すことが多い。 p(x; θ) ​ パラメータが複数ある場合でもまとめてθ と表すことがあるので注意。本資料でも頻繁に出てくる。 10

ガウス分布のパラメータ推定 尤度 観測値xが得られた時のパラメータθ の尤もらしさp(x; θ)。 見ての通り、実態はただの確率密度関数である。ただ尤度を考える場合はパラメータの方を変数と捉える。それを明示するために L(θ) = p(x; θ) ​ と置く場合もある（本資料では使わない） 。 また複数の観測値X = {x(1) , x(2) , ⋯ , x(N ) }に対しては積を取る。 N p(X; θ) = ∏ p(x(n) ; θ) ​ ​ n=1 11

ガウス分布のパラメータ推定 最尤推定 統計では、得られた標本から母集団が従う確率分布のパラメータを推定することがある。標本をよく観察し、適切な分布を仮定し、そのパ ラメータを推定する。 パラメータの推定には様々な手法があり、最も基本的な手法は最尤推定である。尤度が最も大きくなるようなパラメータを推定する。 θ^ = arg max p(X; θ) ​ θ ​ 12

ガウス分布のパラメータ推定 色々なパラメータ推定手法 他にも、事後確率が最大となるようなパラメータを推定するMAP推定や、パラメータも確率変数とみなしてその分布を求めるベイズ推定な どがある。分布を求めるベイズ推定に対し、具体的なパラメータの値を求める最尤推定とMAP推定は点推定と括られることもある。 推定手法 求めるもの 分類 最尤推定 arg maxθ p(X; θ) 点推定 MAP推定 arg maxθ p(θ∣X) 点推定 ベイズ推定 p(θ∣X) ベイズ推定 ​ ​ 13

ガウス分布のパラメータ推定 ガウス分布の最尤推定 ガウス分布のパラメータを最尤推定で推定してみよう。 ある標本X が得られ、以下のようなヒストグラムになったとする。左右対称で中心に多く集中しているため、ガウス分布が仮定できそう。 14

ガウス分布のパラメータ推定 ガウス分布の最尤推定 ガウス分布のパラメータは平均と分散である。直感的には、標本の平均と分散がそのまま入る気がする。実際にそれは正しいが、ここでは 知らないフリをし、ちゃんと計算する。 やり方はシンプルで、尤度関数を微分して傾きが0になる点を求めるだけ。 ただ実際は尤度関数ではなく対数尤度で計算する。その方が計算しやすいから。logは単調増加の関数なので以下が成り立つ。 arg max p(X; θ) = arg max log p(X; θ) ​ θ ​ θ ​ 15

ガウス分布のパラメータ推定 ガウス分布の最尤推定 やってみよう。対数尤度を変形していく。 N log p(X; θ) = log ∏ p(x(n) ; θ) ​ n=1 N = ∑ log p(x(n) ; θ) ​ ​ n=1 N 1 (x(n) − μ)2 = ∑ log ( exp (− )) 2 2 2σ 2πσ n=1 ​ ​ ​ ​ ​ N 1 (x(n) − μ)2 = − ∑ ( log(2π) + log σ + ) 2 2 2σ n=1 ​ ​ ​ これを微分する。 16

ガウス分布のパラメータ推定 ガウス分布の最尤推定 まずμから。対数尤度はμについて上に凸の二次関数となっているた め、唯一の極値が最大値を取る。 N ∂ ∂ (x(n) − μ)2 log p(X; θ) = − ∑ ∂μ ∂μ 2σ 2 n=1 ​ ​ ​ N x(n) − μ =∑ σ2 n=1 ​ ​ これが0になる点を求めると、標本X の平均になる。 ​ ​ ​ N x(n) − μ ∑ =0 σ2 n=1 ​ ​ N ∑ x(n) − N μ = 0 ​ ​ ​ n=1 N 1 ∑ x(n) = μ N n=1 ​ ​ 17

ガウス分布のパラメータ推定 ガウス分布の最尤推定 分散σ 2 でも同じことをやる。微分して N ∂ ∂ ∂ (x(n) − μ)2 log p(X; θ) = − ∑ ( log σ + ) ∂σ ∂σ ∂σ 2σ 2 n=1 ​ ​ ​ ​ ​ N 1 (x(n) − μ)2 = −∑( − ) σ σ3 n=1 ​ ​ N = −∑( ​ n=1 0になる点を求めると、標本X の分散になる。 2 σ について、対数尤度はこの点で唯一の極値をとり、またこの点に おける二階微分が負であるためこの点は極大値である。よってこの点 が最大値である。 ​ ​ ​ σ 2 − (x(n) − μ)2 ) σ3 ​ N σ 2 − (x(n) − μ)2 −∑( )=0 σ3 n=1 ​ ​ N −N σ + ∑(x(n) − μ)2 = 0 2 ​ ​ ​ n=1 N 1 ∑(x(n) − μ)2 = σ 2 N n=1 ​ ​ 18

ガウス分布のパラメータ推定 ガウス分布の最尤推定 ということで、μ ^, σ ^ 2 が推定できた。 ​ N 1 μ ^= ∑ x(n) N n=1 ​ ​ ​ N ​ 1 σ ^2 = ∑(x(n) − μ)2 N n=1 ​ ​ ​ 実際に先の標本の平均と分散を当てはめてみると、良くフィットし ていることが分かる。 19

ガウス分布のパラメータ推定 多変量ガウス分布の最尤推定 1変量と同じ。標本の平均ベクトルと共分散行列になる。 μ ^= ​ 1 ∑ x(n) N n ​ ​ ^ = 1 ∑(x(n) − μ)(x(n) − μ)T Σ N n ​ ​ ​ ​ 導出は省略。 20

混合ガウス分布 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 21

混合ガウス分布 混合ガウス分布 複数のガウス分布を足し合わせた分布。多峰性が見られる分布に適 用する。 K p(x; θ) = ∑ πk Nk (x) ​ ​ ​ ​ k=1 Nk (x) = N (x; μk , σk2 ) πk ∈ [0, 1], ∑k πk = 1 （このπ は円周率じゃないよ） ​ ​ ​ ​ ​ ​ パラメータは、各ガウス分布の平均μk 、分散σk2 と混合比率πk ​ ​ ​ （多変量の場合は共分散行列） 22

混合ガウス分布 混合ガウス分布 K 個のガウス分布N1 , N2 , ⋯ , NK と混合比率π1 , π2 , ⋯ , πK の重み付き和で表されている。 ​ ​ ​ ​ ​ ​ （左がそれぞれのガウス分布で、右がそれを足し合わせたもの） 23

混合ガウス分布 混合ガウス分布を用いたクラスタリング 得られたデータに対し混合ガウス分布を仮定し、パラメータを推定する。その下で以下を計算し、各サンプルx(n) に対しクラスc(n) を割り 当てる。 c(n) = arg max πk Nk (x(n) ) ​ ​ k ​ ​ x(n) について、どのガウス分布から得られる確率が一番高いかを計算する。 24

混合ガウス分布のパラメータ推定 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 25

混合ガウス分布のパラメータ推定 混合ガウス分布の対数尤度 通常のガウス分布では、対数尤度を微分して0になる点を求めるという方法で最尤推定を行なった。混合ガウス分布においても同じやり方 で最尤推定ができないかと考えてみる。対数尤度を見てみよう。 log p(x; θ) = log (∑ πk Nk (x)) K ​ ​ ​ ​ k=1 log-sumと呼ばれるlogの中に∑が入る形になった。この形は非常に厄介で、解析的に解を得られないことが知られている。 こまった、同じ手法が使えない。 26

混合ガウス分布のパラメータ推定 潜在変数 ここで、潜在変数というものを導入してみる。潜在変数とは、観測データの裏に潜む、データをよく表した観測できない変数。 潜在変数の導入とは、観測データX = {x(1) , x(2) , ⋯ , x(N ) }の裏に何らかの潜在変数Z = {z (1) , z (2) , ⋯ , z (N ) }が存在すると仮定 し、それをモデル化すること。xと関係のある何らかの確率変数を考え、それがどのように観測データに作用するかをモデル化するという こと。 27

混合ガウス分布のパラメータ推定 潜在変数 潜在変数の導入を数学的に表現すると、それは観測データと潜在変数の同時分布を定めることになる。 p(x, z) = p(x∣z)p(z) ​ 潜在変数の分布p(z)と、それがデータxにどう作用するかを表す条件付き分布p(x∣z)を定めている。ベイズ的に言うなら事前分布と尤度 関数。 同時分布を周辺化して尤度を表す。 p(x; θ) = ∑ p(x, z; θ) ​ ​ z いきなりこんなことを言われても良くわからないと思うので、例を見てみよう。 28

混合ガウス分布のパラメータ推定 潜在変数 ランダムに日本人を適当な人数選んで身長を計測したとする。そして、この標本から日本人の身長の分布を推定したいとする。ここで潜在 変数を導入してみる。 まず、身長に影響を与えていそうな他の要素を考えてみる。 性別 親の身長 などなど これらは観測できるかもしれないが、身長しか計測していなければ見ることができないので、観測できないという立場をとる。 29

混合ガウス分布のパラメータ推定 潜在変数 次にこの中から潜在変数を仮定する。ここでは性別のみを考慮し、それを潜在変数z ∈ {0, 1}とする。z は男女という2種類のクラスを表 すカテゴリ変数。 ここから同時分布を定義する。まず事前分布だが、z は2値をとる変 数なのでこれはベルヌーイ分布になるね。 次に尤度関数だが、各クラス（性別）に対応するガウス分布としてお こうか。 p(z; θ) = Bernoulli(z; ϕ) p(x∣z; θ) = { ​ N0 (x) N1 (x) ​ ​ ​ ​ z=0 z=1 ​ ​ = Nz (x) ​ これで同時分布を定義できる。 p(x, z; θ) = Bernoulli(z; ϕ)Nz (x) ​ こんな感じで潜在変数を導入する。 30

混合ガウス分布のパラメータ推定 潜在変数 ちなみに、この同時分布を周辺化したものは混合ガウス分布である。峰が二つ見られた分布にK = 2の混合ガウス分布を仮定し、潜在変 数を導入した、という状況。実際、日本人の身長の分布は以下のようになっており、これらの仮定は妥当なものと言える。 引用: 統計情報リサーチ, 男女の身長の分布, https://statresearch.jp/BMI/ 31

• https://statresearch.jp/BMI/

混合ガウス分布のパラメータ推定 混合ガウス分布における潜在変数 混合ガウス分布では、サンプルがどのガウス分布から得られたかを表すカテゴリ変数を潜在変数とする。K 個のガウス分布の中で、サンプ ルx(n) がどこから得られたかを表す潜在変数z (n) ∈ {1, 2, ⋯ , K}を導入する。先の例とほぼ同じで、z が表すクラスが2個からK 個にな っただけ。 p(z; θ) = πz ​ ​ 同時分布は以下。 p(x, z; θ) = p(x∣z; θ)p(z; θ) = Nz (x)πz ​ ​ ​ ​ 実際にこれを周辺化すると先のモデルと一致する。 K K ∑ p(x, z; θ) = ∑ πk Nk (x) ​ z=1 ​ ​ ​ ​ k=1 32

混合ガウス分布のパラメータ推定 潜在変数が与えられた場合 さてここで、仮に潜在変数Z = {z (1) , z (2) , ⋯ , z (N ) }が与えられれば、混合ガウス分布の最尤推定は容易である。各k について、k 番目 のガウス分布から得られたデータXk = {x(n) ∣z (n) = k}を使ってNk (x)の最尤推定を行い、混合比率πk にはNk /N ​ ​ ​ ​ (Nk = ∣Xk ∣)を設 ​ ​ 定すれば良いからね。 しかし、先程も述べた通り、潜在変数Z は観測できない変数であるため、実現できない。 33

混合ガウス分布のパラメータ推定 事後分布が与えられた場合 次に、具体的なz (n) の値ではなく、x(n) が得られた下でのz の確率分布が与えられたと考える。x(n) の事後分布とも見られる。 p(z (n) ) = p(z∣x(n) ) ​ ​ この場合の最尤推定を考える。 34

混合ガウス分布のパラメータ推定 事後分布が与えられた場合 (n) 与えられた確率分布から、観測値x(n) がk 番目のガウス分布から得られた確率が分かる。これを負担率と呼び、γk と置く。 ​ γk(n) = p(z = k ∣x(n) ) ​ p(x(n) , z = k) p(x(n) ) πk Nk (x(n) ) = ∑k πk Nk (x(n) ) = ​ ​ ​ ​ ​ ​ ​ ​ ​ x(i) = 0.5, k = 2の場合の例。 負担率 = (z = k の長さ) / (全体の長さ) 35

混合ガウス分布のパラメータ推定 事後分布が与えられた場合 (n) 負担率γk の和をNk とすると、 ​ ​ Nk = ∑ γk(n) ​ ​ ​ ​ n 以下のようにパラメータが推定できる。 1 (n) ∑ γk x(n) Nk n ^ k = 1 ∑ γ (n) (x(n) − μ )(x(n) − μ )T Σ k k Nk n k μ ^k = ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ π ^k = ​ Nk N ​ ​ (n) 平均、分散をγk の加重平均で推定している。天下り的だが、直感的に受け入れて欲しい。 ​ 36

混合ガウス分布のパラメータ推定 混合ガウス分布のパラメータ推定 さて、こうして得られた3つのパラメータをθ^とし、またそれらによるガウス分布を ^k (x) = N (x; μ ^ k) N ^ k, Σ ​ ​ ​ ​ ​ とすると、新たな負担率を以下に定義できる。 (n) γ^k = p(z = k∣x(n) ; θ^) ^k (x(n) ) π ^k N = ^k (x(n) ) ∑k π ^k N ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ (n) そして、この負担率γ ^k を用いて再度パラメータを推定し…と繰り返すことで対数尤度を大きくすることが可能である。 ​ ​ （証明はEMアルゴリズムの章で。） 37

混合ガウス分布のパラメータ推定 混合ガウス分布のパラメータ推定 ちなみに、クラスタリングで用いるのが負担率である。推定したパラメータで負担率を求め、最も大きな値を取ったk を割り当てる。 (n) c^(n) = arg max γk ​ ​ ​ k 38

混合ガウス分布のパラメータ推定 まとめ 最後にアルゴリズムをまとめる。以下のアルゴリズムで混合ガウス分布の最尤推定ができる。 (n) 1. γk の初期化 ​ 2. θ の推定 (n) 3. γk の更新 ​ 4. 2, 3の繰り返し パラメータθ の初期化から始めてもいい。というか、だいたいそうする。 39

混合ガウス分布のパラメータ推定 おまけ1. k-meansとの関係 (n) 先のアルゴリズムでやっていることはk-meansに近い。γk を0か1に制限した上で重心（平均）μk のみを考慮したものがk-meansである。 ​ ​ c^(n) = arg min ∥x(n) − μk ∥ ​ ​ k 1 γk(n) = { k = c^(n) otherwise 0 1 (n) μ ^k = ∑ γk x(n) Nk n ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ 40

混合ガウス分布のパラメータ推定 おまけ2. EMアルゴリズム 先のアルゴリズムはEMアルゴリズムを混合ガウス分布に適用したものである。EMアルゴリズムとは、混合ガウスモデルのような潜在変数 を仮定できるモデルのパラメータを推定する、より一般的なアルゴリズムである。EMアルゴリズムから先の推定手法を導くこともでき て、それについては最後にまとめた。 先のアルゴリズムは最尤推定を行うものだが、EMアルゴリズムから考えるとMAP推定もできる。また変分ベイズという手法を使うとベイ ズ推定もできる。混合ガウスモデルのパラメータ推定手法をまとめると以下になる 点推定 EMアルゴリズム ベイズ推定 変分ベイズ scikit-learnでもこの二つが用意されている。 sklearn.mixture.GaussianMixture sklearn.mixture.BayesianGaussianMixture 今回は実装するのはEMアルゴリズムによる最尤推定（上のやつ） 。 41

• https://scikit-learn.org/stable/modules/generated/sklearn.mixture.GaussianMixture.html#sklearn.mixture.GaussianMixture
• https://scikit-learn.org/stable/modules/generated/sklearn.mixture.BayesianGaussianMixture.html#sklearn.mixture.BayesianGaussianMixture

混合ガウスモデルの実装 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 42

混合ガウスモデルの実装 混合ガウスモデルの実装 実装してみましょう。必要な数式を次ページ以降にまとめたので、それを元に各自チャレンジしてください。 43

混合ガウスモデルの実装 サンプルデータ 混合ガウスモデルを実装し、これをクラスタリングしよう。 from sklearn.datasets import make_blobs import matplotlib.pyplot as plt d = 2 K = 3 N = 500 X, y = make_blobs( n_samples=N, n_features=d, centers=K, cluster_std=1.5, random_state=9, ) plt.scatter(X[:, 0], X[:, 1]); 44

混合ガウスモデルの実装 テンプレ 以下を完成させよう。 class GaussianMixture: def __init__(self, n_components, max_iter=100): self.n_components = n_components self.max_iter = max_iter def fit(self, X): # 学習部分 def predict(self, X): # 推論部分 # その他必要なメソッド（あれば） n_components がK 。 max_iter が繰り返す回数。 45

混合ガウスモデルの実装 テンプレ こう使って こんな図が出てきたらヨシ gmm = GaussianMixture(n_components=K) gmm.fit(X) clusters = gmm.predict(X) plt.scatter(X[:, 0], X[:, 1], c=clusters); ちなみに、素直に実装しても色が上手く分かれなかったりエラーを吐くことがある。何度か実行してみて、内何度かで図のように3色が綺 麗に分かれればおそらく正解。 46

混合ガウスモデルの実装 実装に必要な情報1. モデル K ∑ πk N (x; μk , Σk ) ​ ​ ​ ​ ​ k=1 K : クラスタ数 パラメータ μ k ∈ Rd ​ Σk ∈ Rd×d πk ∈ [0, 1] ( ∑k πk = 1 ) (k = 1, 2, ⋯ , K) ​ ​ ​ 多変量ガウス分布には scipy.stats.multivariate_normal を使うと良い。 47

• https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.multivariate_normal.html

混合ガウスモデルの実装 実装に必要な情報2. 学習 負担率ではなくパラメータの初期化から始めよう。また、収束判定は考えず、指定した回数の繰り返しとしよう。 1. パラメータの初期化 平均だけはランダムで初期化したほうがいい。 2. 負担率を求める (n) γk = p(z = k∣x(n) , θ) ​ p(x(n) , z = k; θ) = p(x(n) ; θ) πk N (x(n) ; μk , Σk ) = ∑k πk N (x(n) ; μk , Σk ) ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ 48

混合ガウスモデルの実装 実装に必要な情報2. 学習 3. パラメータの更新 μk = ​ 1 (n) ∑ γk x(n) Nk n ​ ​ ​ ​ Σk = ​ ​ 1 (n) ∑ γk (x(n) − μk )(x(n) − μk )T Nk n ​ ​ ​ ​ ​ ​ ​ πk = ​ Nk N ​ ​ ただし Nk = ∑ γk (n) ​ ​ ​ ​ n 4. 2, 3を指定した回数繰り返す 49

混合ガウスモデルの実装 実装に必要な情報3. 推論（クラスタリング） (n) c(n) = arg max γk ​ ​ ​ k 50

混合ガウスモデルの実装 やってみよう！ 次に答え。 51

混合ガウスモデルの実装 実装例 class GaussianMixture: def __init__(self, n_components, max_iter=100): self.n_components = n_components self.max_iter = max_iter def _init_params(self, d): ... def _posterior_dist(self, X): ... def fit(self, X): ... def predict(self, X): ... 1ページに収まらないのでメソッドごとに説明する。全てのコードは当プロジェクトのGitHub参照。 52

混合ガウスモデルの実装 _init_params def _init_params(self, d): self.mu = np.random.randn(self.n_components, d) # (K, d) self.sigma = np.tile(np.eye(d), (self.n_components, 1, 1)) # (K, d, d) self.pi = np.ones(self.n_components) / self.n_components # (K,) パラメータの初期化。K 個をまとめて管理する。 mu sigma pi μ Σ π 形状: (K, d) 形状: (K, d, d) 形状: (K,) ランダムで初期化する。 np.zeros などを用いて全て 全て単位行列で初期化する 一様分布とする 同じ値で初期化してしまうと学習が進まないので注意 53

混合ガウスモデルの実装 _posterior_dist def _posterior_dist(self, X): joint_probs = [] for mu, sigma, pi in zip(self.mu, self.sigma, self.pi): likelihood = sp.stats.multivariate_normal.pdf(X, mean=mu, cov=sigma) joint_prob = pi * likelihood joint_probs.append(joint_prob) joint_probs = np.array(joint_probs) gamma = joint_probs / joint_probs.sum(axis=0) # (K, N) return gamma 事後分布を表す関数。負担率γ を得るための関数。 尤度 likelihood と事前確率 pi の積を同時確率 joint_prob とする。全てのk での同時確率を求め、周辺化したもので割って事後分布 (n) を得る。各x(n) についてK 個のγk をまとめて返すので、このメソッドの出力自体が事後分布になっていると見たほうがいいかもしれな ​ い。 54

混合ガウスモデルの実装 fit def fit(self, X): N, d = X.shape self._init_params(d) # 1. パラメータ初期化 for _ in range(self.max_iter): # 4. 繰り返し gamma = self._posterior_dist(X) # 2. 負担率の計算 Nk = gamma.sum(axis=1) # 3. パラメータ（mu）の更新 self.mu = gamma @ X / Nk.reshape(-1, 1) # 3. パラメータ（Sigma）の更新 dev = np.array([X - self.mu[k] for k in range(self.n_components)]) dev_gamma = (gamma[..., np.newaxis] * dev).transpose(0, 2, 1) self.sigma = dev_gamma @ dev / Nk.reshape(-1, 1, 1) # 3. パラメータ（pi）の更新 self.pi = Nk / N 積の和は行列積を使うと上手く求められるので、μとΣの更新ではそれを使って和を求め、最後にNk で割っている。 ​ 55

混合ガウスモデルの実装 predict def predict(self, X): gamma = self._posterior_prod(X) c = gamma.argmax(axis=0) return c 負担率が最も高くなったクラスk を出力。 56

混合ガウスモデルの実装 上手くいかないケース1 先の実装例は理論を再現できているが、完璧なものではなく、学習が上手くいかないケースが存在する。 一つ目は、計算誤差により共分散行列Σが正則にならないケース。何度か実行すると稀に発生する。逆行列が存在しないため、多変量ガウ ス分布で計算ができず、 scipy.stats.multivariate_normal でエラーが出る。 57

混合ガウスモデルの実装 上手くいかないケース2 二つ目は、クラスタリングが上手くいかないケース。 後の章でも述べるが、先で実装したアルゴリズムは局所解に陥るといった状況が起こり得る。そうなると上手くいかない。これはしょうが ない。初期値次第。 scikit-learnでは初期値をk-meanによって求めることでこのケースが少なくなるよう工夫している。 58

EMアルゴリズム 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 59

EMアルゴリズム EMアルゴリズム 潜在変数を仮定したモデルの点推定を行う数値的なアルゴリズム。混合ガウスモデルの学習アルゴリズムよりも一般的なもの。 ここでは、EMアルゴリズムによる最尤推定を見ていく。 60

EMアルゴリズム EMアルゴリズム 最尤推定なので、対数尤度の最大化という目標は変わらない。ということで対数尤度を見ていくのだが、ここで、イェンゼンの不等式とい うものがある。 イェンゼンの不等式 g(x)が上に凸の関数であるとき、任意のxi と、λi ≥ 0, ∑i λi = 1を満たす任意のλi で以下が成り立つ。 ​ ​ ​ ​ ​ ∑ λi g(xi ) ≤ g (∑ λi xi ) ​ i ​ ​ ​ ​ ​ ​ i 61

EMアルゴリズム EMアルゴリズム イェンゼンの不等式を使って対数尤度を変形する。 log p(X; θ) = log ∑ p(X, Z; θ) ​ Z = log ∑ q(Z) ​ Z ​ p(X, Z; θ) q(Z) 1 ← q(z) q(z) = をかけただけ ​ ​ p(X, Z; θ) ≥ ∑ q(Z) log q(Z) Z ​ ​ ​ = L(q(Z), θ) 途中の不等式がイェンゼンの不等式から得られたもの。q(Z)はイェンゼンの不等式を使うために適当に持ち出したなんらかの確率分布。 現時点では明確になっていない。 最後に出てきた式は下界と呼び、L(q(Z), θ)と置いた。エビデンス（対数尤度）の下界ということで、ELBO（Evidence lower bound）と も呼ぶ。変分ベイズでも出てくるので、変分下界（Variational Lower Bound）とも呼ぶ。 62

EMアルゴリズム EMアルゴリズム ここで、対数尤度と下界の差を求めると log p(X; θ) − L(q(Z), θ) p(X, Z; θ) q(Z) = log p(X ; θ) − ∑ q (Z ) log ​ ∑z q(z) = 1より、 log p(X; θ) = ∑Z q(Z) log p(X; θ) ​ ​ Z → = ∑ q(Z) log p(X; θ) − ∑ q(Z) log ​ ​ Z Z ​ p(X∣Z; θ)p(Z; θ) q(Z) ​ = ∑ q(Z)( log p(X; θ) − log p(X∣Z; θ) − log p(Z; θ) + log q(Z)) ​ ​ ​ Z ベイズの定理 → = ∑ q(Z)( − log ​ Z p(X∣Z; θ)p(Z; θ) + log q(Z)) p(X; θ) ​ = ∑ q(Z)( − log p(Z∣X; θ) + log q(Z)) ​ Z = DKL [q(Z)∣∣p(Z∣X; θ)] ​ 63

EMアルゴリズム EMアルゴリズム q(Z)とp(Z∣X; θ)のKLダイバージェンスになった。KLダイバージェンスは二つの確率分布の近さを表す指標として知られる。 これを以下のようにまとめる。 log p(X; θ) = L(q(Z), θ) + DKL [q(Z)∣∣p(Z∣X; θ)] ​ ​ 対数尤度をKLダイバージェンスと下界の和で表すことが出来た。 64

EMアルゴリズム EMアルゴリズム 我々はこれを最大化したい。ここで思いつくのは、下界とKLダイバージェンスを別々に最大化する方法である。そうすれば対数尤度も大き くなりそう。しかし、このやり方は上手くいかない。 我々が自由に操作できる変数は以下の二つである。 θ q(Z) この二つの変数を対数尤度が最大になるように設定したい訳だが、よく見ると、対数尤度はq(Z)に依存していない。つまり、q(Z)を変化 させた時に変わるものは対数尤度ではなく、下界とKLダイバージェンスの比率である。どちらかを大きくするq(Z)は、もう一方を小さく するq(Z)となる。 65

EMアルゴリズム EMアルゴリズム ではどうすべきか。正解は以下。 1. θ を初期化する 2. θ を固定し、KLダイバージェンスを最小化する 3. q(Z)を固定し、下界を最大化する 4. 2, 3 を繰り返す まずはθ の初期化だがこれは何でもよい。とりあえずθ old と表記しておく。 次に、そのθ old の下でKLダイバージェンスを最小化する。つま りDKL [q(Z)∣∣p(Z∣X; θ old )] = 0となるq(Z)を求める。 ここで、KLダイバージェンスは二つの確率分布が完全に一致する際に唯一の0を ​ とる。ということで、q(Z) = p(Z∣X; θ old )とする。 66

EMアルゴリズム EMアルゴリズム 次にそのq(Z)を固定して下界を最大化する。 θnew = arg max L(p(Z∣X; θold ), θ) ​ ​ θ ここで下限を少し整理する。 L(p(Z∣X; θold ), θ) = ∑ p(Z∣X; θold ) log ​ Z p(X, Z; θ) p(Z∣X; θold ) ​ = ∑ p(Z∣X; θold ) log p(X, Z; θ) + const ​ ​ ​ Z = Ep(Z∣X;θold ) [ log p(X, Z; θ)] + const ​ = Q(θ∣θold ) + const 事後分布に関する同時分布の対数尤度の期待値になった。これはQ(θ, θold )と表しておく。この期待値を最大にするθ を求め、それを新た なθ とする。この期待値はsum-logの形になっており、解析的に解きやすい。 θnew = arg max Q(θ, θold ) ​ θ ​ 67

EMアルゴリズム EMアルゴリズム そしてθ old := θ new として、同じように期待値を最大化する。これを繰り返すことで、対数尤度を大きくすることができる。 ということで、EMアルゴリズムの流れは以下。 1. 潜在変数の導入 2. パラメータθ old の初期化 3. 同時分布の期待値の最大化: θnew ← arg maxθ Q(θ∣θold ) θold ← θnew ​ 4. 3を繰り返す 期待値（Expected value）の最大化（Maximization）を繰り返すため、EM アルゴリズムと呼ばれる。また、数学的には期待値を求めるス テップと最大化のステップに分けられ、それぞれEステップ・Mステップと呼ばれたりする。 68

EMアルゴリズム 対数尤度が単調増加することの証明 以下を示す。 log p(X; θnew ) ≥ log p(X; θold ) 69

EMアルゴリズム 対数尤度が単調増加することの証明 q(Z) = p(Z∣X; θ)であるとき、対数尤度と下界は一致する。よって log p(X; θold ) = L(p(Z∣X; θold ), θold ) ​ ここで、 θnew = arg max L(p(Z ∣X ; θold ), θ) ​ θ ​ より、 L(p(Z∣X; θold ), θnew ) ≥ L(p(Z∣X; θold ), θold ) ​ が成り立つ。 70

EMアルゴリズム 対数尤度が単調増加することの証明 また、任意のq(Z)で log p(X; θ) ≥ L(p(Z), θ) ​ が成り立つので、 log p(X ; θnew ) ≥ L(p(Z ∣X ; θold ), θnew ) ​ が成り立つ。よって log p(X; θnew ) ≥ L(p(Z∣X; θold ), θnew ) ​ ≥ L(p(Z∣X; θold ), θold ) ​ = log p(X; θold ) が成り立つ。証明終わり。 71

EMアルゴリズム 対数尤度が単調増加することの証明 なお、今で示したのは対数尤度が単調増加することで、パラメータが最適解に収束することではない。いわゆる局所解に陥るといった状況 はEMアルゴリズムでも起こり得るので、初期値依存の手法ではある。 72

EMアルゴリズム EMアルゴリズムの可視化 q(Z) = p(Z∣X; θold )としたとき、下界と対数尤度は以下のような関係で、 log p(X; θ) ≥ L(p(Z∣X; θold ), θ) ​ ​ θ = θold の場合にのみイコールとなる。 対数尤度は複雑で解析的に解が得られない。一方で下界はシンプルであるため解析的に解を求められる。ここまでの関係を可視化する。 73

EMアルゴリズム EMアルゴリズムの可視化 EMアルゴリズムではこの下界を最大化し、θ new とする。 74

EMアルゴリズム EMアルゴリズムの可視化 そしてθ old := θ new と更新した上で、同じように下界を最大化する。 これを繰り返すことでパラメータを推定する。 75

EMアルゴリズム EMアルゴリズムの可視化 EMアルゴリズムが初期値に依って最適解に収束しない可能性があることは図からも分かり、この例だと初期値が右側になった場合、右の 峰にハマってしまいそう。 76

EMアルゴリズムの混合ガウスモデルへの適用 1. ガウス分布 2. ガウス分布のパラメータ推定 3. 混合ガウス分布 4. 混合ガウス分布のパラメータ推定 5. 混合ガウスモデルの実装 6. EMアルゴリズム 7. EMアルゴリズムの混合ガウスモデルへの適用 77

EMアルゴリズムの混合ガウスモデルへの適用 EMアルゴリズムの混合ガウス分布への適用 EMアルゴリズムから混合ガウス分布のパラメータ推定手法を導出する。 モデルを再掲。 p(z) = πz ​ p(x∣z) = Nz (x) ​ p(x; θ) = ∑ p(x, z; θ) ​ ​ z K ​ = ∑ πk Nk (x) ​ ​ ​ k=1 78

EMアルゴリズムの混合ガウスモデルへの適用 同時分布の期待値 EMアルゴリズムでは同時分布の期待値の最大化を繰り返す。期待値を変形。 Q(θ, θold ) = Ep(Z∣X;θold ) [ log p(X, Z; θ) ] ​ = Ep(Z∣X;θold ) [∑ log p(x(i) , z (i) ; θ)] ​ ​ i = ∑ Ep(z∣x(i) ;θold ) [ log p(x(i) , z; θ) ] ​ ​ i = ∑ Ep(z∣x(i) ;θold ) [ log Nz (x(i) ) + log πz ] ​ ​ ​ ​ ​ i ​ = ∑ ∑ p(z = k∣x(i) ; θold ) (log Nk (x(i) ) + log πk ) ​ i ​ ​ ​ k = ∑ ∑ γk (log N (x(i) ; μk , Σk ) + log πk ) (i) ​ i ​ ​ ​ ​ ​ k 1 (i) − μk )) + log πk ) = ∑ ∑ γk(i) (− (d log(2π) + log ∣Σ∣k + (x(i) − μk )T Σ−1 k (x 2 i ​ ​ ​ ​ ​ ​ ​ ​ ​ k 79

EMアルゴリズムの混合ガウスモデルへの適用 負担率 これはさっきと一緒。 (i) γk = p(z = k∣x(i) , θold ) ​ p(x(i) ∣z = k; θold )p(z = k; θold ) p(x(i) ; θold ) old πkold N (x(i) ; μold k , Σk ) = old ∑k πkold N (x(i) ; μold k , Σk ) = ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ Nk = ∑ γk(n) ​ ​ ​ ​ ​ n 80

EMアルゴリズムの混合ガウスモデルへの適用 パラメータの最大化: μk ​ ここから、各パラメータで偏微分して0になる点を求めていく。まずμk から。 ​ ∂Q 1 (n) ∂ = − ∑ γk (x(n) − μk )T Σ−1 (x(n) − μk ) ∂μk 2 n ∂μk ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ (n) = −Σ−1 − μk ) k ∑ γk (x (n) ​ ​ ​ ​ ​ n 81

EMアルゴリズムの混合ガウスモデルへの適用 パラメータの最大化: μk ​ 0になる点を求める (n) −Σ−1 − μk ) = 0 k ∑ γk (x (i) ​ ​ ​ ​ n 1 (n) ∑ γk x(n) = μk Nk n ​ ​ ​ ​ ​ ​ ​ 82

EMアルゴリズムの混合ガウスモデルへの適用 パラメータの最大化: Σk 次、Σk ​ ​ 偏微分 ∂Q 1 (n) ∂ (n) ( log ∣Σk ∣ + (x(n) − μk )T Σ−1 = − ∑ γk − μk )) k (x ∂Σk 2 n ∂Σk 1 (n) −1 (n) = − ∑ γk (Σ−1 − μk )(x(n) − μk )T Σ−1 k − Σk (x k ) 2 n ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ 0になる点を求める 83

EMアルゴリズムの混合ガウスモデルへの適用 パラメータの最大化: πk ​ πk は他のパラメータと違い、制約条件がある。 ​ g(πk ) = ∑ πk′ − 1 = 0 ​ ​ ​ ​ k′ ラグランジュの未定乗数法を使おう。 ∂(Q + λg(πk )) 1 (n) ∂ ( log πk + λ( ∑ πk′ − 1)) = − ∑ γk ∂πk 2 n ∂πk k′ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ 1 = ∑ γk(n) ( + λ) πk n ​ ​ ​ ​ ​ 84

EMアルゴリズムの混合ガウスモデルへの適用 パラメータの最大化: πk 0になる点を求める。 ​ ∑ γk ( (n) ​ ​ n 1 + λ) = 0 πk ​ ​ − ∑ γk = λπk (n) ​ ​ ​ ​ ​ n ここで、両辺でk について総和を取り、 − ∑ ∑ γ k = ∑ λ πk (n) ​ ​ ​ n k ​ ​ k ​ ​ −N = λ これを代入する。 − ∑ γk = −N πk (n) ​ ​ ​ n 1 (n) ∑ γ = πk N n k ​ ​ ​ ​ ​ ​ 85

EMアルゴリズムの混合ガウスモデルへの適用 まとめ μk = ​ 1 ∑ γk(n) x(n) Nk n ​ ​ ​ ​ Σk = ​ ​ 1 (n) ∑ γk (x(n) − μk )(x(n) − μk )T Nk n ​ ​ ​ ​ ​ ​ ​ πk = ​ Nk N ​ ​ さっきの解が得られた。 86

EMアルゴリズムの混合ガウスモデルへの適用 参考 1. 後藤 正幸, 小林 学, 入門 パターン認識と機械学習, コロナ社 2. 石塚 崚斗, Academaid, 【徹底解説】EMアルゴリズムをはじめからていねいに, https://academ-aid.com/ml/em, 参照: 2024-05-31 3. 斎藤 康毅, ゼロから作るDeep Learning ❺ ―生成モデル編, O’Reilly Japan 87

• https://academ-aid.com/ml/em

オワリ おつ