這一節學習內容為概率圖模型裡的一節，因為下午在跑程式手裡也沒什麼事情幹，寫個em的demo記錄一下。本文也不是來推導過程的，只是方法和**記錄，推導請看其他部落格。
概率圖模型跳不過的一章就是高斯混合模型，高斯混合模型是由多個高斯分布以一定權重組成的模型，其概率密度函式等於各高斯模型的概率密度函式加權求和。
一般地，每個高斯分布的概率密度為：
\[n(x|\mu_k,\sigma_k) = \frac\sigma} exp(-\frac)
\]其中k代表第k個高斯分布，因此gmm的概率分布為：
\[p(x) = \sum_^k \pi_kn(x|\mu_k,\sigma_k)
\]其中\(\sum_^ = 1\)且\(\pi_k \geq 0\).
gmm的乙個基本問題就是選用高斯混合模型對資料進行建模，那麼模型的引數如何獲得呢？
乙個gmm的引數僅由一組\(、、(\pi、\mu、\sigma)\)決定，所以我們要估計的引數就是這三個。其涉及的問題是包含隱變數的引數的求解過程，如果模型不含隱變數，直接用最大似然然後求個導就能解決，但是現在模型含有隱變數，也就意味著求導的過程極其複雜，乙個比較好的方法是通過em演算法來進行引數估計。
關於gmm的簡介就到這裡，網上有很多關於gmm的介紹，也非常通俗易懂，不多贅述。
em演算法是乙個求解帶隱變數的模型的引數常用的方法，其通過最大化訓練集的邊界似然，迭代更新引數。
過程為：
先初始化隨機引數
e步：固定引數求解後驗概率
m步：固定後驗概率，優化求解證據下界對應的引數。
重複進行直到收斂。
對於gmm模型而言，要求其引數就是按照上面的步驟進行。但在開始之前，我們要先生成給定分布的樣本，即取樣。取樣步驟為：
按照\(、\pi_1、\pi_2 \dots \pi_k\)的的分布隨機選擇乙個高斯分布
假設選擇了第k個高斯分布，按照\(n(x|\mu_k,\sigma_k)\)的分布取樣乙個樣本(高斯分布的取樣計算機可以實現)
固定\(\mu,\sigma\)，計算後驗概率分布\(p(z^|x^)\)
\[\gamma_ \xlongequal p(z^ = k|x^) \\=\frac)p(x^|z^)})} \\=\frac)p(x^|z^)}^\pi_kn(x^|\mu_k,\sigma_k)}
\]用上式\(\gamma_\)表示第n個樣本對第k個高斯分布的後驗概率，其用\(n \times k\)的矩陣表示。
在已知後驗概率的情況下，最大化邊際似然\(p(x|\mu,\sigma)\),即最大化對數似然:
\[log(p(x|\mu,\sigma)) = \sum_^ log \sum_^ q(z)\frac,z^=k|\mu_k,\sigma_k)} = k)}\\\geq \sum_^\sum_^ log \space q(z)\frac,z^=k|\mu_k,\sigma_k)} = k)}\\\xlongequal elbo(q,x|\mu,\sigma) \\
\]其中\(q(z) = p(z^=k)\),elbo通過不等式定義為函式下界，所以我們要求的其實是對數似然的下界。
\[\max_ elbo(\gamma,d|\pi,\mu,\sigma)\\s.t. \sum_^ = 1
\]d為訓練集，\(\gamma\)為後驗分布。
利用拉格朗日法求解後有如下更新的結論：
\[\pi_k \leftarrow \frac\\\mu_k \leftarrow \frac\sum_^n\gamma_x^\\\sigma_k = \sqrt\sum_^n \gamma_(x^)^2}
\]其中\(n_k = \sum_^n \gamma_\)
'''
@descripttion: this is aoru xue's demo, which is only for reference.
@version:
@author: aoru xue
@date: 2019-12-30 20:46:28
@lasteditors : aoru xue
@lastedittime : 2019-12-31 00:03:01
'''from matplotlib import pyplot as plt
import numpy as np
import random
import pickle
# **裡的迴圈運算可以矩陣優化，有興趣且需要請自行優化，高維需要修改部分**。為保證例子簡單，故以一維高斯分布為例。
class gmm():
def __init__(self,pi,mu,sigma,samples): # 生成乙個一維高斯混合分布
self.x = np.zeros(shape = (samples,))
for i in range(samples):
idx = np.random.choice(range(pi.shape[0]),p=pi)
self.x[i] = np.random.normal(mu[idx],sigma[idx])
self.mu = np.random.randn(mu.shape[0])
self.sigma = np.ones(sigma.shape[0])
self.pi = np.ones(sigma.shape) / sigma.shape[0]
@staticmethod
def gausian_probablity(x,mu,sigma):
return 1. / np.sqrt(2 * np.pi)/sigma * np.exp(- (x - mu)**2 / 2 / sigma**2)
def em_method(self,):
post_probability = np.zeros((self.x.shape[0],self.pi.shape[0]),dtype = np.float) # (n,k)
for _ in range(100):
# efor n in range(post_probability.shape[0]):
for k in range(post_probability.shape[1]):
post_probability[n,k] = self.gausian_probablity(self.x[n],self.mu[k],self.sigma[k]) * self.pi[k]
post_probability = post_probability / np.sum(post_probability,axis = 1,keepdims=true) # 後驗概率
EM演算法學習筆記詳解
基於李航的 統計學習方法 對em演算法學習的介紹性筆記。一.預備知識 在正式介紹em演算法之前，先介紹推導em演算法用到的數學基礎知識，包括凹凸函式，jensen不等式。上凸函式 函式f x 滿足對定義域上任意兩個數a,b都有f a b 2 f a f b 2 下凸函式 函式f x 滿足對定義域上任...
EM演算法學習筆記 1 對EM演算法的簡單理解
因做實驗的需要，最近在學習 em演算法，演算法介紹的資料網上是有不少，可是沒有一篇深入淺出的介紹，演算法公式太多，比較難懂，畢竟她是 ml領域 10大經典演算法之一 且一般是結合 gmm模型的引數估計來介紹 em的。看過不少 em的資料，現將自己對 em演算法用稍微通俗點的文字寫下來，當然你可以用 ...
機器學習演算法 EM演算法
介紹 理論分析 現在,已經知道目標函式 似然函式 maxj p x 假設隱變數為z,那麼可以寫成j zp x z 我們假設當前的引數為 i 那麼可以得到 j j i logp x logp x i log zp x,z logp x i log z p z x,i p x z p z x,i log...