目标检测 PAA - 高斯混合模型（GMM）和期望最大化算法（EM algorithm）

目标检测 PAA - 高斯混合模型（GMM）和期望最大化算法（EM algorithm）

flyfish

论文 Probabilistic Anchor Assignment with IoU Prediction for Object Detection

关键字:
期望最大化算法(EM算法) expectation maximization algorithm
高斯混合模型( GMM) Gaussian Mixture Model

我们想拥有这样的数据集，有身高数据项也有性别数据项

但是实际情况我们有这样的数据集，只有身高，但是没有性别。

问题来了
在一个学校里共有500多个学生，我们选择200个学生，100个男生的身高和100个女生的身高，共200条记录，但是我们不知道这200个数据中哪个是男生的身高，哪个是女生的身高。假设男生、女生的身高分别服从正态分布，不知道每个样本即每条记录从哪个分布抽取的。对于每一个样本，就有两个方面需要估计： 这个身高数据是关于男生还是女生？男生、女生身高的正态分布的参数分别是多少？这个问题可以直接跳到下面场景3问题的解决

有一个数据集 一半数据来自高斯分布A，一半数据来自高斯分布B，两种数据的来源都是不同的高斯分布随机生成的，能够描述这个数据集的模型就是高斯混合模型（Gaussian Mixed Model）
这个模型不是来自一个分布所以是混合的,所以就叫了高斯混合模型。

我们根据此提出三个问题
场景1
数据未标注，已知参数。未标注的数据表示我们不知道一个样本是来自分布A，还是分布A，但是我们知道两个高斯分布的参数即期望和标准差（mean、standard deviation 没有上标2,所以不是方差）

场景2
数据已标注，未知的参数。我们知道每个样本是来自分布A是分布B，但是分布 AB的参数即期望和标准差我们都不知道

场景3
数据未标注，参数也未知。这种场景怎么处理呢？

场景1问题的解决

第一个场景用概率密度函数（probability density function）来解决问题

公式可以不知道，只需要知道两个参数 $\mu ,\sigma$
一元高斯高分布（univariate Gaussian distribution）
$$p\left(x\mid \mu ,{\sigma }^2\right)=\mathcal{N}\left(\mu ,{\sigma }^2\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}\exp \left(-\frac{(x-\mu {)}^2}{2{\sigma }^2}\right).$$

多元高斯分布（二元）（multivariate Gaussian distribution ）
$$p(\mathbf{x}\mid \mathbf{\mu },\mathbf{\Sigma })=\mathcal{N}\left(\mathbf{\mu },\mathbf{\Sigma }\right)=(2\pi {)}^{-\frac{D}{2}}|\mathbf{\Sigma }{|}^{-\frac{1}{2}}\exp \left(-\frac{1}{2}(\mathbf{x}-\mathbf{\mu }{)}^{\top }{\mathbf{\Sigma }}^{-1}(\mathbf{x}-\mathbf{\mu })\right).$$
根据两个参数我们可视化下,将值带入公式就行


实现代码

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import multivariate_normal
from mpl_toolkits.mplot3d import Axes3D
"""
返回高斯分布数据
参数:
x -- 输入
mu -- 均值
sigma -- 标准差
# 根据公式变代码
"""
def Gaussian_Data(x, mu, sigma):
    t1=2*(sigma**2)
    exp = np.exp( -(np.square(x - mu) / t1) )
    t2 = 1 / (np.sqrt(2 * np.pi) * sigma)
    r = t2 * exp
    return r

#一元高斯高分布
mu=0.04780295283011753
sigma=1.1251987967998018
x = np.linspace(-1,1,100)
y = Gaussian_Data(x,mu,sigma)
plt.plot(x, y)
plt.show()

#多元高斯分布（二元）
# 生成二维网格平面
X, Y = np.meshgrid(np.linspace(-1,1,100), np.linspace(-1,1,100))
# 二维坐标数据
d = np.dstack([X,Y])
#标量变矩阵
mu = np.zeros(2) + mu  # 均值矩阵
sigma = np.eye(2) * sigma  # 协方差矩阵

Z = multivariate_normal(mean=mu, cov=sigma).pdf(d)#用库函数
fig = plt.figure(figsize=(6,4))
ax = Axes3D(fig)
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='seismic', alpha=0.8)
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

场景2问题的解决

第二个场景用 最大似然估计（Maximum Likelihood estimation）来解决问题
场景2的代码实现

'''
norm.cdf 返回对应的累计分布函数值
norm.pdf 返回对应的概率密度函数值
norm.rvs 产生指定参数的随机变量
norm.fit 返回给定数据下，各参数的最大似然估计（MLE）值
'''
from scipy.stats import norm
import matplotlib.pyplot as plt
import numpy as np

x_norm = norm.rvs(size=100)
#在这组数据下，正态分布参数的最大似然估计值
x_mean, x_std = norm.fit(x_norm)#Maximum Likelihood Estimate.
print ('mean, ', x_mean)
print ('x_std, ', x_std)
plt.hist(x_norm,density=True)
x = np.linspace(-5,5,100)
plt.plot(x, norm.pdf(x), 'r-')#Probability density function evaluated at x
plt.show()

假设学校所有学生的身高服从高斯分布 $N(\mu ,{\sigma }^2)$ 。不知道分布的参数即均值 $\mu$ 和方差 ${\sigma }^2$ ，如果我们想估计出这两个参数，那么问题就解决了。那么怎样估计这两个参数呢？
假设我们随机选择 200 个人。然后统计随机抽取的这 200 个人的身高，根据这 200 个人的身高估计均值 $\mu$ 和方差 ${\sigma }^2$。
为了计算学校学生的身高分布，概率密度 $p(x|\theta )$ 服从高斯分布 $N(\mu ,{\sigma }^2)$ ,从中随机选择 200 个样本，组成样本集 X = { x 1 , x 2 , … , x N } X=\{x_1,x_2,…,x_N\} X={x1​,x2​,…,xN​} ，其中 $x_i$ 表示选择的第 $i$ 个人的身高，N 等于 200，表示样本个数，我们想通过样本集 $X$来估计出未知参数 $\theta$。这里概率密度 $p(x|\theta )$ 服从高斯分布 $N(\mu ,{\sigma }^2)$ ，其中的未知参数是 $\theta =(\mu ,\sigma )$。
怎样估算参数 $\theta$ 呢？

由于每个样本都是独立地从 $p(x|\theta )$ 中抽取的，他们之间没有关系所以相互独立的。假设选择学生 A的概率是 $p(x_A|\theta )$，选择学生B的概率是 $p(x_B|\theta )$，那么同时抽到男生 A 和男生 B 的概率是 $p(x_A|\theta )\times p(x_B|\theta )$，同时随机选择这 200 个学生的概率就是他们各自概率的乘积，即为联合概率 。
$\prod$这个符号表示连乘
$$L(\theta )=L(x_1,x_2,\cdots ,x_n;\theta )=\prod _{i=1}^{n}p(x_i|\theta ),\theta \in \Theta$$

函数反映的是在不同的参数 $\theta$ 取值下，取得当前这个样本集的可能性，因此称为参数 $\theta$ 相对于样本集 X 的似然函数（likelihood function），记为 $L(\theta )$。

为什么是这200 个身高数据，而不是其他，因为这 200 个身高数据出现的概率极大，也就是其对应的似然函数 $L(\theta )$ 极大。
就像一个骰子有6个面，其中有5个面都是6点，那么仍骰子大概率是6点。
$$\hat{\theta }=\text{argmax}L(\theta )$$
$\hat{\theta }$ 这个叫做 $\theta$ 的极大似然估计，即为我们所求的值。
利用最大似然估计计算出均值和标准差.
到这里是不是感觉到了交叉熵。

场景3问题的解决

也就是知道数据标注可以求出参数，知道参数也可以求出数据标注，都是一个已知，求另一个未知，两个都未知怎么求呢？也就是我们开头提出的学生问题，俗称鸡和蛋的问题。
我们用实际数据来演示此问题是如何解决的
先将数据可视化下，在已知男女标签和身高的情况下的可视化
数据集从我的github上下载，实际我们不知道gender那列都是什么数据。


可视化代码

import numpy as np
import seaborn as sns
from matplotlib import pyplot as plt

data = np.genfromtxt('./Expectation_Maximization.csv', delimiter=',', skip_header=1)
print(data.shape)
x = data[0:-1, -2]
N=247
x1 = x[:N]
x2 = x[N:]
#可视化1
bins = np.linspace(150, 200, 100)
plt.hist(x1, bins, alpha=0.5, label='x1',density=True)
plt.hist(x2, bins, alpha=0.5, label='x2',density=True)
plt.legend(loc='upper right')
plt.show()

#可视化2
fig, ax = plt.subplots()
ax.set_xlim([150, 200])
sns.distplot(x1, color='green')
sns.distplot(x2, color='orange')
plt.show()

以上两个场景相当于已知x计算y和已知y计算x
我们虽然不知道，但是我们可以随机初始化两个高斯分布。如果随机选择的参数恰好是我们要求的值，这个小概率事件，所以还得有下步
我们可以使用随机参数初始化两个高斯分布，然后在两个步骤之间进行迭代：
（一）估计标签，同时保持参数固定（第一个场景），
（二）更新参数，同时保持标签固定（第二个场景），
在这两个步骤上进行迭代最终将达到局部最优。

一个固定，更新另一个 像不像多元函数求导数的方法即求偏导数
还使用迭代方式像极了梯度下降。
上述的迭代算法就是EM 算法，全称 Expectation Maximization Algorithm。期望最大算法是一种迭代算法，用于含有隐变量（Hidden Variable）的概率参数模型的最大似然估计。
代码如下

import numpy as np
from scipy.stats import norm
import scipy.stats as stats
import matplotlib.pyplot as plt
import matplotlib


def plot_distributions(data, data_sampled, mu, sigma, K, color="green", color_sampled="red", name='plot.png'):
    matplotlib.rcParams['text.usetex'] = True
    plt.rcParams.update({'font.size': 16})
    data_sampled = np.clip(data_sampled, np.min(data), np.max(data))
    plt.hist(data, bins=15, color=color, alpha=0.45, density=True)
    plt.hist(data_sampled, bins=15, range=(np.min(data), np.max(data)), color=color_sampled, alpha=0.45, density=True)
    for k in range(K):
        curve = np.linspace(mu[k] - 10 * sigma[k], mu[k] + 10 * sigma[k], 100)
        color = np.random.rand(3)
        plt.plot(curve, stats.norm.pdf(curve, mu[k], sigma[k]), color=color, linestyle="--", linewidth=3)
    plt.ylabel(r"$p(x)$")
    plt.xlabel(r"$x$")
    plt.tight_layout()
    #plt.xlim(20, 120)
    plt.xlim(150, 200)
    plt.savefig(name, dpi=200)
    plt.show()


def plot_likelihood(nll_list):
    matplotlib.rcParams['text.usetex'] = True
    plt.rcParams.update({'font.size': 16})
    plt.plot(np.arange(len(nll_list)), nll_list, color="black", linestyle="--", linewidth=3)
    plt.ylabel(r"(negative) log-likelihood")
    plt.xlabel(r"iteration")
    plt.tight_layout()
    plt.xlim(0, len(nll_list))
    plt.savefig('nll.png', dpi=200)
    plt.show()


def sampler(pi, mu, sigma, N):
    data = list()
    for n in range(N):
        k = np.random.choice(len(pi), p=pi)
        sample = np.random.normal(loc=mu[k], scale=sigma[k])
        data.append(sample)
    return data


def main():
    data = np.genfromtxt('./Expectation_Maximization.csv', delimiter=',', skip_header=1)  # [:,-2]
    print(data.shape)
    data = data[0:500, -2]
    print(data)
    # print(data.shape,data[0:5 , -1])
    # return
    N = data.shape[0]
    K = 2  # two components GMM
    tot_iterations = 100  # stopping criteria

    # Step-1 (Init)
    mu = np.random.uniform(low=160.0, high=200.0, size=K)
    sigma = np.random.uniform(low=20.0, high=30.0, size=K)
    # mu = np.random.uniform(low=42.0, high=95.0, size=K)
    # sigma = np.random.uniform(low=5.0, high=10.0, size=K)
    pi = np.ones(K) * (1.0 / K)  # mixing coefficients
    r = np.zeros([K, N])  # responsibilities
    nll_list = list()  # store the neg log-likelihood

    for iteration in range(tot_iterations):
        # Step-2 (E-Step)
        for k in range(K):
            r[k, :] = pi[k] * norm.pdf(x=data, loc=mu[k], scale=sigma[k])
        r = r / np.sum(r, axis=0)  # [K,N] -> [N]

        # Step-3 (M-Step)
        N_k = np.sum(r, axis=1)  # [K,N] -> [K]
        for k in range(K):
            # update means
            mu[k] = np.sum(r[k, :] * data) / N_k[k]
            # update variances
            numerator = r[k] * (data - mu[k]) ** 2
            sigma[k] = np.sqrt(np.sum(numerator) / N_k[k])
            # update weights
        pi = N_k / N

        likelihood = 0.0
        for k in range(K):
            likelihood += pi[k] * norm.pdf(x=data, loc=mu[k], scale=sigma[k])
        nll_list.append(-np.sum(np.log(likelihood)))
        # Check for invalid negative log-likelihood (NLL)
        # The NLL is invalid if NLL_t-1 < NLL_t
        # Note that this can happen for round-off errors.
        if (len(nll_list) >= 2):
            if (nll_list[-2] < nll_list[-1]): raise Exception("[ERROR] invalid NLL: " + str(nll_list[-2:]))

        print("Iteration: " + str(iteration) + "; NLL: " + str(nll_list[-1]))
        print("Mean " + str(mu) + "\nStd " + str(sigma) + "\nWeights " + str(pi) + "\n")

        # Step-4 (Check)
        if (iteration == tot_iterations - 1): break  # check iteration

    plot_likelihood(nll_list)
    data_gmm = sampler(pi, mu, sigma, N=100)
    plot_distributions(data, data_gmm, mu, sigma, K, color="green", color_sampled="red", name="plot_sampler.png")


if __name__ == "__main__":
    main()

我这里没做的工作是将500多条数据shuffle,然后再选择200条数据。
现在sklearn代码库已经包含GaussianMixture对象，
GaussianMixture 对象实现了用来拟合高斯混合模型的 期望最大化 (EM) 算法。它还可以为多变量模型绘制置信椭圆体，同时计算 BIC（Bayesian Information Criterion，贝叶斯信息准则）来评估数据中聚类的数量。而其中的GaussianMixture.fit 方法可以从训练数据中拟合出一个高斯混合模型。
如果给定测试数据，通过使用 GaussianMixture.predict 方法，可以为每个样本分配最适合的高斯分布模型。库的使用方法如下

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.colors import LogNorm
from sklearn import mixture
from mpl_toolkits.mplot3d import Axes3D
n_samples = 300

# generate random sample, two components
np.random.seed(0)

# generate spherical data centered on (20, 20)
shifted_gaussian = np.random.randn(n_samples, 2) + np.array([20, 20])

# generate zero centered stretched Gaussian data
C = np.array([[0., -0.7], [3.5, .7]])
stretched_gaussian = np.dot(np.random.randn(n_samples, 2), C)

# concatenate the two datasets into the final training set
X_train = np.vstack([shifted_gaussian, stretched_gaussian])

# fit a Gaussian Mixture Model with two components
clf = mixture.GaussianMixture(n_components=2, covariance_type='full')
clf.fit(X_train)

# display predicted scores by the model as a contour plot
x = np.linspace(-20., 30.)
y = np.linspace(-20., 40.)
X, Y = np.meshgrid(x, y)
XX = np.array([X.ravel(), Y.ravel()]).T
Z = -clf.score_samples(XX)
Z = Z.reshape(X.shape)

CS = plt.contour(X, Y, Z, norm=LogNorm(vmin=1.0, vmax=1000.0),
                 levels=np.logspace(0, 3, 10))
#CB = plt.colorbar(CS, shrink=0.8, extend='both')
plt.scatter(X_train[:, 0], X_train[:, 1], .8)
plt.clabel(CS,fontsize=10,colors=('r'),fmt='%.2f')
plt.title('Negative log-likelihood predicted by a GMM')
plt.axis('tight')
plt.show()

问题解了，解的就是高斯混合模型 Gaussian Mixture Model (GMM)

EM算法可以参考
《机器学习与应用》雷明
16章介绍 高斯混合模型 P474
18章介绍 EM算法即期望最大化算法 P508
其他参考
What is the expectation maximization algorithm? Chuong B Do & Serafim Batzoglou
mixture models
categorical distribution
probability distribution
Kalman filter
Gaussian processes
determinant
likelihood
independent and identically distributed
law of large numbers
Density Estimation for a Gaussian mixture