机器学习代码实现2：贝叶斯模型

这一专题的内容主要是我在进行统计机器学习的算法的代码实践的过程中的一些记录，包括公式的推导，算法的实现与数据集的实验结果。

这部分工作我比较早之前就已经开始做了，现在刚好趁此机会整理一下之前写过的所有小demo，汇总成一个完成的库

1.贝叶斯定理

  • 贝叶斯定理主要是用来判断一个数据样本是否属于类别，其理论主要就是基于下面的公式：

  • 这里的称为先验概率，也就是目前已经知道的类别的分布情况，而叫做似然likelihood，表示的是当前已知的属于类别的的情况，也就是说我们要判断一个样本是否属于类别，我们需要先得到已经属于类别的同类的概率和类别本身出现的概率，这样就可以估计当前的属于类别的概率
  • 当有C个类别的时候，贝叶斯定理可以表示为：

  • 计算出结果之后，根据极大似然的规则，让可能性最高的作为样本类别的预测结果。
  • 因此在训练一个贝叶斯分类器的时候，我们要得到一个类别-特征的分布矩阵，我们可以假设可能的预测结果有类，而不同的特征可能取值有种，这样一来，我们就可以得到一个的矩阵，矩阵中的每一个元素代表了某一种特征情况输入这一类别的概率(或者出现的次数也可以，需要再进一步处理)

2.贝叶斯定理代码实现

  • 代码实现中我们重点需要关注和的计算，对于，我们可以对矩阵先按行求和，得到训练集中每种样本的出现数量，然后除以样本总数就可以得到每种类别对应的概率分布，然后是对于矩阵的每一行中的元素分别乘以该行的元素总量之和
  • 最终的贝叶斯定理代码实现如下所示：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
class Bayes:
  	# 计算似然，也就是P(x_i|\omega_j)的值
    def likelihood(self, X: np.ndarray) -> np.ndarray:
        """
        计算数据集的似然
        :param X: 输入的数据集，表示不同类别的特征分布情况，大小是C*N，其中C是类别的个数，N是特征的个数
        :return: 数据集的似然
        """
        self.C, self.N = X.shape
        self.likelihood = np.zeros(X.shape)
        class_sum = np.sum(X, axis=1)
        for i in range(self.C):
            self.likelihood[i, :] = X[i, :] / class_sum[i]
        return self.likelihood
     # 计算后验概率，当然这里写的比较naive没有向量化
    def posterior(self, X: np.ndarray) -> np.ndarray:
        self.C, self.N = X.shape
        self.p = np.zeros(X.shape)
        likelihood = self.likelihood(X)
        prior = np.sum(X, axis=1) / np.sum(X)
        p_X = np.zeros(self.N)
        for j in range(self.N):
            for i in range(self.C):
                p_X[j] += likelihood[i, j] * prior[i]
        for i in range(self.C):
            for j in range(self.N):
                self.p[i, j] = likelihood[i, j] * prior[i] / p_X[j]
        return self.p

3.数据集实验

  • 这里使用的数据集来自浙江大学蔡登教授开设的《机器学习》课程，这里的数据集仅有1维特征和标签，比较简单。我们分别对仅使用似然和使用后验概率来进行预测的结果进行了可视化，得到的结果是：

    • 似然：

    似然分布情况

    • 后验概率：

    后验概率分布情况

4.参数估计与代码实现

  • 事实上我们也发现了贝叶斯定理的核心问题在于如何计算，这里我们可以使用正态分布作为工具，多远正态分布的定义式是：

  • 而正态分布模型的参数估计的表达式是：

  • 我们可以用样本数据集X来估计出正态分布模型的参数，然后用这个正态分布模型来计算likelihood，再根据贝叶斯定理算出后验概率，最后按照最大的概率作为决策。
  • 参数估计的代码实现如下：
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
def gaussian_pos_prob(X, Mu, Sigma, Phi):
    """
    GAUSSIAN_POS_PROB Posterior probability of GDA.
    Compute the posterior probability of given N data points X
    using Gaussian Discriminant Analysis where the K gaussian distributions
    are specified by Mu, Sigma and Phi.
    Inputs:
        'X'     - M-by-N numpy array, N data points of dimension M.
        'Mu'    - M-by-K numpy array, mean of K Gaussian distributions.
        'Sigma' - M-by-M-by-K  numpy array (yes, a 3D matrix), variance matrix of
                  K Gaussian distributions.
        'Phi'   - 1-by-K  numpy array, prior of K Gaussian distributions.
    Outputs:
        'p'     - N-by-K  numpy array, posterior probability of N data points
                with in K Gaussian distribsubplots_adjustutions.
    """
    N = X.shape[1]
    K = Phi.shape[0]
    p = np.zeros((N, K))
    # 先计算likelihood
    likelihood = np.zeros((N, K))
    for i in range(N):
        p_x = 0
        for j in range(K):
            x_minus_mu = X[:, i] - Mu[:, j]
            sigma = Sigma[:, :, j]
            # 估计参数，并构造一个高斯函数
            det_sigma = np.linalg.det(sigma)
            inv_sigma = np.linalg.inv(sigma)
            base = 1.0 / (2 * np.pi * np.sqrt(np.abs(det_sigma)))
            exponent = np.matmul(np.matmul(x_minus_mu.T, inv_sigma), x_minus_mu) * -0.5
            # 计算似然
            likelihood[i, j] = base * np.exp(exponent)
            p_x += likelihood[i, j] * Phi[j]
        for j in range(K):
            p[i, j] = likelihood[i, j] * Phi[j] / p_x
    return p