[人工智能]（七）高斯判别分析

前言

这篇文章总结一个生成模型——高斯判别分析。

1. 判别模型和生成模型

?在监督学习中，模型训练成功后，我们往往会得到一个决策函数 $f(x)$ 或者一个概率分布 $P(y|x)$ ，当需要预测时，我们通过给定的 $x$ ，使用 $f(x)$ 作为 $y$ 的预测值或者使用 $P(y|x)$ 计算出概率。

?在之前总结的线性回归、逻辑回归等模型都是判别模型(discriminative model)，判别模型是根据训练集学习，直接得到 $f(x)$ 或者 $P(y|x)$ ，它只是关心对于给定的 $x$ 需要预测出什么样的 $y$ 。

?还有一种模型是生成模型(generative model)，这个模型大致思想是这样的，对于给定的训练集 ${\left\{\left(x^{\left(i\right)},y^{\left(i\right)}\right)\right\}}_1^m$ ，我们它是由某一个概率分布 “生成” 的，这样我们就要估计出那个概率分布。这些数据既然会出现，那么也就是说出现这些数据的可能性很大，我们就通过最大化联合概率分布 $P\left(x^{\left(1\right)},y^{\left(1\right)},?,x^{\left(m\right)},y^{\left(m\right)}\right)$ 来估计原来的分布。也就是利用最大似然估计法估计。估计原分布的过程就是训练的过程。预测时，通过条件概率公式 $P\left(y|x\right)=P\left(x,y\right)/P\left(x\right)$ 得到给定特征向量 $x$ 条件下 $y$ 的概率。

2. 高斯判别分析模型

?高斯判别分析(Gaussian Discriminant Analysis，GDA) 是一种生成模型。假设一个分类任务中，类别有 $K$ 个，记为 { 0 , 1 , ? ? , K ? 1 } \left\{ 0,1,\cdots ,K-1 \right\} {0,1,?,K?1} ，特征向量为 $x={\left(x_1,x_2,?,x_n\right)}^T$ ，输出变量为 y ∈ { 0 , 1 , ? ? , K ? 1 } y\in \left\{ 0,1,\cdots ,K-1 \right\} y∈{0,1,?,K?1} ，现有训练集 ${\left\{\left(x^{\left(i\right)},y^{\left(i\right)}\right)\right\}}_1^m$ 。

?高斯判别分析需对每一个类别 $i$ 作出这样的假设：
$$P\left(x|y=i\right)=\frac{1}{{\left(2\pi \right)}^{n/2}{\left|\Sigma \right|}^{1/2}}\exp \left(?\frac{1}{2}{\left(x?{\mu }_i\right)}^T{\Sigma }^{?1}\left(x?{\mu }_i\right)\right)$$
也就是说，每一类别 $i\left(i=0,1,?,K?1\right)$ 的特征向量服从 $n$ 元高斯分布。${\mu }_i$ 是 $n$ 维期望向量，第 $j$ 个分量是类别 $i$ 的特征向量 $x$ 的第 $j$ 个分量 $x_j$ 的期望。 $\Sigma$ 是 $n\times n$ 维协方差矩阵，${\Sigma }_{ij}$ 是 $x_i$ 和 $x_j$ 的协方差，$\Sigma$ 并不区分类别，这里假设所有类别的 $\Sigma$ 相同。

?需要对类别的概率分布作出这样的假设：
P ( y ) = ? 0 I { y = 0 } ? 1 I { y = 1 } ? ? K ? 1 I { y = K ? 1 } P\left( y \right) =\phi _{0}^{I\left\{ y=0 \right\}}\phi _{1}^{I\left\{ y=1 \right\}}\cdots \phi _{K-1}^{I\left\{ y=K-1 \right\}} P(y)=?0I{y=0}??1I{y=1}???K?1I{y=K?1}?
其中 ${\sum }_{i=0}^{K?1}{?}_i=1$ ，也就是说 $y$ 是类别 $i\left(i=0,1,?,K?1\right)$ 的概率为 ${?}_i$ 。

?这个模型训练时就是要求解这些参数： { ? 0 , ? 1 , ? ? , ? K ? 1 , μ 0 , μ 1 , ? ? , μ K ? 1 , Σ } \left\{ \phi _0,\phi _1,\cdots ,\phi _{K-1},\mu _0,\mu _1,\cdots ,\mu _{K-1},\Sigma \right\} {?0?,?1?,?,?K?1?,μ0?,μ1?,?,μK?1?,Σ} ，实际上，对于每一个类别很可能有不同的协方差矩阵 $\Sigma$ ，但这个模型中仍然假设 $\Sigma$ 相同，这样做可以使得分类决策面是线性的。但是如果每个类别的协方差矩阵相差很大时，那么这样做可能使模型预测效果变差。

?模型预测时，对于所给的特征向量 $x$ ，计算出条件概率：
$$\begin{gathered} P\left(y=i|x\right)=\frac{P\left(x,y=i\right)}{P\left(x\right)}=\frac{P\left(x|y=i\right)P\left(y=i\right)}{P\left(x\right)} \\ \left(i=0,1,?,K?1\right) \end{gathered}$$
?预测结果是 $y=arg{\max }_{i}P\left(y=i|x\right)$ 。

?无论 $i$ 取值如何，概率 $P(x)$ 都是相同的，因此预测时可以不必计算 $P(x)$ ，只需要计算 $P(x|y=i)P(y=i)$ ，然后求出哪一类 $i$ 使得 $P(x|y=i)P(y=i)$ 最大即可。这也就是说预测结果可以简化为：$y=arg{\max }_{i}P\left(y=i|x\right)=arg{\max }_{i}P\left(x|y=i\right)P\left(y=i\right)$ 。

?如果想要得到每个类别的概率，这时只需要将每一个 $P(x|y=i)P(y=i)$ 归一化即可，也就是说：
$$P\left(y=i|x\right)=\frac{P\left(x|y=i\right)P\left(y=i\right)}{{\sum }_{j=0}^{K?1}P\left(x|y=j\right)P\left(y=j\right)}\left(i=0,1,?,K?1\right)$$
?实际上，通过全概率公式也可以得到：$P\left(x\right)={\sum }_{j=0}^{K?1}P\left(x|y=j\right)P\left(y=j\right)$ 。

3. 参数估计

?模型中提到，训练时要估计这些参数： { ? 0 , ? 1 , ? ? , ? K ? 1 , μ 0 , μ 1 , ? ? , μ K ? 1 , Σ } \left\{ \phi _0,\phi _1,\cdots ,\phi _{K-1},\mu _0,\mu _1,\cdots ,\mu _{K-1},\Sigma \right\} {?0?,?1?,?,?K?1?,μ0?,μ1?,?,μK?1?,Σ} ，采用最大似然法估计，对数似然函数为：
$$\begin{gathered} l\left({\mu }_0,?,{\mu }_{K?1},{?}_0,?,{?}_{K?1},\Sigma \right) \\ =\sum _{i=1}^m\ln P\left(x^{\left(i\right)},y^{\left(i\right)}\right)=\sum _{i=1}^m\ln P\left(x^{\left(i\right)}|y^{\left(i\right)}\right)P\left(y^{\left(i\right)}\right) \end{gathered}$$
最大化该对数自然函数，求出的结果为：
$${\mu }_i=\frac{{\sum }_{j=1}^m{x}^{\left(j\right)}I\left\{y^{\left(j\right)}=i\right\}}{{\sum }_{j=1}^{m}I\left\{y^{\left(j\right)}=i\right\}}\left(i=0,1,?,K?1\right)$$

$${?}_i=\frac{1}{m}\sum _{j=1}^{m}I\left\{y^{\left(j\right)}=i\right\}$$

$$\Sigma =\frac{1}{m}\sum _{j=1}^m\left(x^{\left(j\right)}?{\mu }_{y^{\left(j\right)}}\right){\left(x^{\left(j\right)}?{\mu }_{y^{\left(j\right)}}\right)}^T$$

?最大似然估计最终得出的结果是上述公式。这个计算过程是稍微有些繁琐的（所以我也就不想写了），但是这个最终结果却是非常棒的，看起来很合理，直接使用矩估计也可以得出同样的结果。

4. 代码实现

?实现高斯判别分析的类：

# GDA类实现
class GDA(object):
    # GDA的参数
    mu, Sigma, phi, n_class = 0, 0, 0, 0

    def fit(self, X_train, y_train):
        self.n_class = y_train.max() + 1 # 总共有这么多类别，从 0 编号
        # 计算 mu
        self.mu = np.empty([self.n_class, X_train.shape[1]])
        for i in range(self.n_class):
            idx = np.argwhere(y == i)
            idx = idx.squeeze(axis = 1)
            self.mu[i] = X[idx].mean(axis = 0)
        # 计算协方差矩阵
        self.Sigma = np.corrcoef(X_train, rowvar = False)
        # 计算输出变量分布的参数
        self.phi = np.empty(self.n_class)
        for i in range(self.n_class):
            self.phi[i] = (y_train == i).sum() / y_train.shape[0]
    
    def probability(self, X):
        mu, Sigma, phi, n_class = self.mu, self.Sigma, self.phi, self.n_class
        p = np.empty([X.shape[0], n_class]) # 求解概率
        for i in range(X.shape[0]):
            for j in range(n_class):
                # 分子部分（多元高斯分布密度函数的分子）
                fenzi = np.exp(-np.dot(np.dot(X[i]-mu[j], np.linalg.inv(Sigma)), X[i]-mu[j]) / 2) 
                # 下面是密度函数的分母部分（这部分可以不要，因为最终我们是会归一化的，分母只是个常数）
                # fenmu = (2 * np.pi)**(X.shape[1] / 2) * np.linalg.det(Sigma)**(1 / 2) 
                p[i,j] = phi[j] * fenzi
        p = (p.T / p.sum(axis = 1)).T
        return p # 返回这个概率
    
    def predict(self, X):
        mu, Sigma, phi, n_class = self.mu, self.Sigma, self.phi, self.n_class
        p = np.empty([X.shape[0], n_class]) # 求解概率的主要部分（不包括常数部分）
        for i in range(X.shape[0]):
            for j in range(n_class):
                tmp = np.exp(-np.dot(np.dot(X[i]-mu[j], np.linalg.inv(Sigma)), X[i]-mu[j]) / 2) # 分子部分
                p[i, j] = phi[j] * tmp
        pred = np.argmax(p, axis = 1) # 最大概率就是预测结果
        return pred # 返回预测结果

? 自己创造数据并使用高斯判别分析训练并预测：

if __name__ == '__main__':
    # 创造数据
    # 这个是二分类数据
    X, y = make_blobs(n_samples = 200, n_features = 2, centers = 2, cluster_std = (1,1), random_state = 0)
    # 这个是三分类数据
    # X, y = make_blobs(n_samples = 200, n_features = 2, centers = 2, cluster_std = (1,1), random_state = 0)
    
    plt.scatter(X[:,0], X[:,1], c = y, s = 6) # 可视化
    X_train, X_test, y_train, y_test = train_test_split(
        X, y, train_size = 0.7, random_state = 0) # 训练集，测试集
    # 测试模型
    clf = GDA()
    clf.fit(X_train, y_train)
    y_pred = clf.predict(X_test) # 预测
    y_prob = clf.probability(X_test) # 概率
    conf_mat = confusion_matrix(y_test, y_pred) # 混淆矩阵
    print("混淆矩阵为：\n", conf_mat)
    
    # 画出区域的预测图
    x1, x2 = np.mgrid[X[:,0].min() - 1:X[:,0].max() + 1:0.05,X[:,1].min() - 1:X[:,1].max() + 1:0.05]
    area = np.stack([x1.flat,x2.flat],axis=1)
    area_pred = clf.predict(area)
    area_pred = area_pred.reshape(x1.shape)
    plt.pcolormesh(x1, x2, area_pred, alpha = 0.2)
    plt.show()

• 创造的二分类数据为：

模型划分测试集与训练集，训练后在测试集上预测，预测结果的混淆矩阵为：

模型在区域上的分类结果为：

• 创造的三分类数据为：

模型训练后在测试集上预测，预测结果的混淆矩阵为：

模型在区域上的预测效果为：

?从上述两个实例中可以看到，模型大致成功分类了样本，且分类决策面是线性的，这是因为 $\Sigma$ 选取一致。

创作打卡挑战赛

赢取流量/现金/CSDN周边激励大奖