贝叶斯算法

贝叶斯定理相关公式

先验概率P(A)：在不考虑任何情况下，A事件发生的概率。

条件概率P(B|A)：A事件发生的情况下，B事件发生的概率。

后验概率P(A|B)：在B事件发生之后，对A事件发生的概率的重新评估。

全概率：如果A和A'构成样本空间的一个划分，那么事件B的概率为：A和A'的概率分别乘以B对这两个事件的概率之和。

基于条件概率的贝叶斯定律数学公式

朴素贝叶斯算法

朴素贝叶斯(Naive Bayes， NB)是基于“特征之间是独立的”这一朴素假设，应 用贝叶斯定理的监督学习算法

对应给定的样本X的特征向量x1 ,x2 ,...,xm；该样本X的类别y的概率可以由贝叶斯公式得到：

特征属性之间是独立的，所以得到：

公式优化得到：

在给定样本的情况下，P(x₁ ,x₂ ,...,x_m)是常数，所以得到：

从而：

朴素贝叶斯算法流程

朴素贝叶斯算法流程/定义如下：

• 设x={x1 ,x2 ,...,xm}为待分类项，其中xi为x的一个特征属性。
• 类别集合为C={y1 ,y2 ,...,yn }。
• 分别计算P(y1 |x),P(y2 |x),....,P(yn |x)的值（贝叶斯公式）。
• 如果P(yk |x)=max{P(y1 |x),P(y2 |x),....,P(yn |x)},那么认为x为yk类型。

高斯朴素贝叶斯

Gaussian Naive Bayes是指当特征属性为连续值时，而且分布服从高斯分布，那么在计算P(x|y)的时候可以直接使用高斯分布的概率公式：

因此只需要计算出各个类别中此特征项划分的各个均值和标准差。

伯努利朴素贝叶斯

Bernoulli Naive Bayes是指当特征属性为连续值时，而且分布服从伯努利分布，那么在计算P(x|y)的时候可以直接使用伯努利分布的概率公式：

伯努利分布是一种离散分布，只有两种可能的结果。1表示成功，出现的概率为p； 0表示失败，出现的概率为q=1-p；其中均值为E(x)=p，方差为Var(X)=p(1-p)。

多项式朴素贝叶斯

Multinomial Naive Bayes是指当特征属性服从多项分布(特征是离散的形式的时候)，直接计算类别数目的占比作为先验概率和条件概率。

N是总样本个数，k是总的类别个数，N_yk是类别为y_k的样本个数，α为平滑值。

N_yk是类别为y_k的样本个数，n为样本的维度数目，N_yk,xi为类别y_k中第i维特征的值为x_i的样本个数，α为平滑值。

当α=1时，称为Laplace平滑，当0<α<1时，称为Lidstone平滑，α=0时不做平滑；平滑的主要作用是可以克服条件概率为0的问题。

多项式朴素贝叶斯案例理解

对于下列训练数据，使用多项式朴素贝叶斯方式对测试样本(2,M)和(3,S)做一个 预测判断。

import numpy as np
import pandas as pd
# x1 = np.array([1, 1, 1, 1, 2, 2, 2, 3, 3, 4])
# x2 = np.array(['S', 'M', 'S', 'L', 'S', 'S', 'L', 'L', 'L', 'S'])
# y = np.array([-1, 1, 1, 1, -1, -1, 1, 1, -1, 1])

x1 = np.array([1, 1, 1, 1, 2, 2, 2, 2, 3, 4])
x2 = np.array(['S', 'M', 'S', 'L', 'S', 'S', 'L', 'L', 'L', 'S'])
y = np.array([-1, 1, 1, 1, -1, -1, 1, -1, 1, 1])
train_D = np.vstack((x1,x2)).T
class MNB():
    def __init__(self,alpha=0):
        self.alpha = alpha
    def fit(self,X, Y):
        X_size = X.shape
        N = X.shape[0]
        n = X_size[1]
        N_yk = {}
        xi_list = []
        for i in range(n):
            xi_list.append(X[:, i])
        self.cls = np.unique(Y.reshape(-1))
        k = len(self.cls)
        for c in self.cls:
            N_yk['N_y=%s' % c] = len(X[y == c])

        N_yk_xi = {}
        for c in self.cls:
            for i, d in enumerate(xi_list):
                for e in np.unique(d):
                    N_yk_xi['Ny=%s,x%s=%s' % (c, i + 1, e)] = len(train_D[(y == c) & (d == e)])

        # 先验概率
        self.P_yk = {}
        for c in self.cls:
            self.P_yk['P_y=%s' % c] = (N_yk['N_y=%s' % c] + self.alpha) / (N + k * self.alpha)

        self.P_yk_xi = {}
        # 条件概率
        for c in self.cls:
            for i, d in enumerate(xi_list):
                for e in np.unique(d):
                    self.P_yk_xi['P_y=%s,x%s=%s' % (c, i + 1, e)] = (N_yk_xi['Ny=%s,x%s=%s' % (c, i + 1, e)] + self.alpha) / (
                                N_yk['N_y=%s' % c] + n * self.alpha)

    def predict(self,datas):
        result = []
        self.predict_proba = []
        for x in datas:
            pre = []
            for c in self.cls:
                p_c = self.P_yk['P_y=%s' % c]
                for i, e in enumerate(x):
                    p_c *= self.P_yk_xi['P_y=%s,x%s=%s' % (c, i + 1, e)]
                pre.append(p_c)
            self.predict_proba.append(pre)
            result.append(self.cls[pre.index(max(pre))])
        self.predict_proba = np.array(self.predict_proba)
        return result

test = np.array([[2,'M'],[3,'S']])
alpha = 1
algo = MNB(alpha)
algo.fit(train_D,y)
print('alpha=%s的预测结果：%s'%(alpha,algo.predict(test)))
print('各个类别的预测概率')
print(algo.predict_proba)


'''基于数组实现概率矩阵'''
class MultiNB():
    import numpy as np
    import pandas as pd
    def __init__(self,alpha=1):
        self.alpha = alpha

    def fit(self,X,Y):
        # class LenError(Exception): #定义一个异常类
        #     pass
        # if len(X) != len(Y):
        #     raise LenError('X与Y的长度不一致')

        assert len(X) == len(Y),'X与Y的长度不一致' #使用断言语句
        X = np.array(X)
        Y = np.array(Y).reshape(-1)
        x_size = X.shape
        N = x_size[0]
        self.n = x_size[1]
        self.Y_k = {i:k for i,k in enumerate(np.unique(Y))}
        self.k = len(self.Y_k)
        NB_P = []
        for c in self.Y_k.values():
            P_dict = {}
            c_datas = X[Y==c]
            N_yk = len(c_datas)
            P_dict['P(yk)'] = (N_yk+self.alpha)/(N+self.k*self.alpha)
            for i in range(self.n):
                x_i = c_datas[:,i]
                for x in np.unique(X[:,i]):
                    N_yk_xi = np.sum(x_i == x)
                    P_dict['P(x%s=%s|yk)'%(i+1,x)] = (N_yk_xi+self.alpha)/(N_yk+self.n*self.alpha)
            NB_P.append(pd.Series(P_dict))
        self.NB_P = pd.DataFrame(NB_P)
        self.NB_P.replace(np.nan,0,inplace=True)
        self.NB_P.index.name = 'yk'
        print(self.NB_P)

    def predict(self,X):
        result = []
        self.predict_proba = []
        for x in X:
            pre = []
            for c in range(self.k):
                P = self.NB_P.loc[c,'P(yk)']
                for i, e in enumerate(x):
                    P *= self.NB_P.loc[c,'P(x%s=%s|yk)'%(i+1,e)]
                pre.append(P)
            self.predict_proba.append(pre)
            result.append(self.Y_k[pre.index(max(pre))])
        result = np.array(result)
        try:
            result = result.astype(np.float)
        except:
            pass
        self.predict_proba = np.array(self.predict_proba)
        return result

#显示被省略的列
pd.set_option('display.width',None)
alog = MultiNB(alpha)
alog.fit(train_D,y)
print('alpha=%s的预测结果：%s'%(alpha,alog.predict(test)))
print('各个类别的预测概率')
print(alog.predict_proba)

sklearn中朴素贝叶斯的案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl
from sklearn.preprocessing import StandardScaler, MinMaxScaler, PolynomialFeatures
from sklearn.naive_bayes import GaussianNB, MultinomialNB#高斯贝叶斯和多项式朴素贝叶斯
from sklearn.pipeline import Pipeline
from sklearn.metrics import accuracy_score
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier

## 设置属性防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False

# 花萼长度、花萼宽度，花瓣长度，花瓣宽度
iris_feature_E = 'sepal length', 'sepal width', 'petal length', 'petal width'
iris_feature_C = [u'花萼长度', u'花萼宽度', u'花瓣长度', u'花瓣宽度','class']
iris_class = 'Iris-setosa', 'Iris-versicolor', 'Iris-virginica'
features = [2,3]

## 读取数据
path = '../datas/iris.data'  # 数据文件路径
data = pd.read_csv(path, header=None)
x = data[list(range(4))]
x = x[features]
y = pd.Categorical(data[4]).codes ## 直接将数据特征转换为0，1,2
print ("总样本数目：%d；特征属性数目：%d" % x.shape)

## 0. 数据分割，形成模型训练数据和测试数据
x_train1, x_test1, y_train1, y_test1 = train_test_split(x, y, train_size=0.8, random_state=14)
x_train, x_test, y_train, y_test = x_train1, x_test1, y_train1, y_test1
print ("训练数据集样本数目：%d, 测试数据集样本数目：%d" % (x_train.shape[0], x_test.shape[0]))

## 高斯贝叶斯模型构建
clf = Pipeline([
        # ('sc', StandardScaler()),#标准化，把它转化成了高斯分布, 特征取值存在负数
        ('poly', PolynomialFeatures(degree=4)),
        ('clf', MultinomialNB())]) # MultinomialNB多项式贝叶斯算法中要求特征属性的取值不能为负数
## 训练模型
clf.fit(x_train, y_train)

## 计算预测值并计算准确率
y_train_hat = clf.predict(x_train)
print ('训练集准确度: %.2f%%' % (100 * accuracy_score(y_train, y_train_hat)))
y_test_hat = clf.predict(x_test)
print ('测试集准确度：%.2f%%' % (100 * accuracy_score(y_test, y_test_hat)))


## 产生区域图
N, M = 500, 500     # 横纵各采样多少个值
x1_min1, x2_min1 = x_train.min()
x1_max1, x2_max1 = x_train.max()
x1_min2, x2_min2 = x_test.min()
x1_max2, x2_max2 = x_test.max()
x1_min = np.min((x1_min1, x1_min2))
x1_max = np.max((x1_max1, x1_max2))
x2_min = np.min((x2_min1, x2_min2))
x2_max = np.max((x2_max1, x2_max2))

t1 = np.linspace(x1_min, x1_max, N)
t2 = np.linspace(x2_min, x2_max, N)
x1, x2 = np.meshgrid(t1, t2)  # 生成网格采样点
print(x1.shape)
x_show = np.dstack((x1.flat, x2.flat))[0] # 测试点

cm_light = mpl.colors.ListedColormap(['#77E0A0', '#FF8080', '#A0A0FF'])
cm_dark = mpl.colors.ListedColormap(['g', 'r', 'b'])
y_show_hat = clf.predict(x_show)                  # 预测值
y_show_hat = y_show_hat.reshape(x1.shape)

## 画图
plt.figure(facecolor='w')
plt.pcolormesh(x1, x2, y_show_hat, cmap=cm_light)     # 预测值的显示
plt.scatter(x_train[features[0]], x_train[features[1]], c=y_train, edgecolors='k', s=50, cmap=cm_dark)
plt.scatter(x_test[features[0]], x_test[features[1]], c=y_test, marker='^', edgecolors='k', s=120, cmap=cm_dark)
plt.xlabel(iris_feature_C[features[0]], fontsize=13)
plt.ylabel(iris_feature_C[features[1]], fontsize=13)
plt.xlim(x1_min, x1_max)
plt.ylim(x2_min, x2_max)
plt.title(u'GaussianNB对鸢尾花数据的分类结果, 正确率:%.3f%%' % (100 * accuracy_score(y_test, y_test_hat)), fontsize=18)
plt.grid(True)
plt.show()

import numpy as np
import pandas as pd
from time import time
import matplotlib.pyplot as plt
import matplotlib as mpl

from sklearn.datasets import fetch_20newsgroups#引入新闻数据包
from sklearn.feature_extraction.text import TfidfVectorizer#做tfidf编码
from sklearn.feature_selection import SelectKBest, chi2#卡方检验——特征筛选
from sklearn.linear_model import RidgeClassifier
from sklearn.svm import LinearSVC,SVC
from sklearn.naive_bayes import MultinomialNB, BernoulliNB, GaussianNB #引入多项式和伯努利的贝叶斯
from sklearn.neighbors import KNeighborsClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import GridSearchCV
from sklearn import metrics

###  基准模型方法
def benchmark(clf,name):
    print (u'分类器：', clf)
    
    ##  设置最优参数，并使用5折交叉验证获取最优参数值
    alpha_can = np.logspace(-2, 1, 10)
    model = GridSearchCV(clf, param_grid={'alpha': alpha_can}, cv=5)
    m = alpha_can.size
    
    ## 如果模型有一个参数是alpha，进行设置
    if hasattr(clf, 'alpha'):
        model.set_params(param_grid={'alpha': alpha_can})
        m = alpha_can.size
    ## 如果模型有一个k近邻的参数，进行设置
    if hasattr(clf, 'n_neighbors'):
        neighbors_can = np.arange(1, 15)
        model.set_params(param_grid={'n_neighbors': neighbors_can})
        m = neighbors_can.size
    ## LinearSVC最优参数配置
    if hasattr(clf, 'C'):
        C_can = np.logspace(1, 3, 3)
        model.set_params(param_grid={'C':C_can})
        m = C_can.size
    ## SVM最优参数设置
    if hasattr(clf, 'C') & hasattr(clf, 'gamma'):
        C_can = np.logspace(1, 3, 3)
        gamma_can = np.logspace(-3, 0, 3)
        model.set_params(param_grid={'C':C_can, 'gamma':gamma_can})
        m = C_can.size * gamma_can.size
    ## 设置深度相关参数，决策树
    if hasattr(clf, 'max_depth'):
        max_depth_can = np.arange(4, 10)
        model.set_params(param_grid={'max_depth': max_depth_can})
        m = max_depth_can.size
    
    ## 模型训练
    t_start = time()
    model.fit(x_train, y_train)
    t_end = time()
    t_train = (t_end - t_start) / (5*m)
    print (u'5折交叉验证的训练时间为：%.3f秒/(5*%d)=%.3f秒' % ((t_end - t_start), m, t_train))
    print (u'最优超参数为：', model.best_params_)
    
    ## 模型预测
    t_start = time()
    y_hat = model.predict(x_test)
    t_end = time()
    t_test = t_end - t_start
    print (u'测试时间：%.3f秒' % t_test)
    
    ## 模型效果评估
    train_acc = metrics.accuracy_score(y_train, model.predict(x_train))
    test_acc = metrics.accuracy_score(y_test, y_hat)
    print (u'训练集准确率：%.2f%%' % (100 * train_acc))
    print (u'测试集准确率：%.2f%%' % (100 * test_acc))
    
    ## 返回结果(训练时间耗时，预测数据耗时，训练数据错误率，测试数据错误率, 名称)
    return t_train, t_test, 1-train_acc, 1-test_acc, name


## 设置属性防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'SimHei']
mpl.rcParams['axes.unicode_minus'] = False


## 数据加载
print (u'加载数据...')
t_start = time()
## 不要头部信息
remove = ('headers', 'footers', 'quotes')
## 只要这四类数据（可以先加载所有数据，然后通过target_names获取取值）
categories = 'alt.atheism', 'talk.religion.misc', 'comp.graphics', 'sci.space'
# categories = None

## 分别加载训练数据和测试数据
data_train = fetch_20newsgroups(data_home='./datas/',subset='train', categories=categories, shuffle=True, random_state=0, remove=remove)
data_test = fetch_20newsgroups(data_home='./datas/',subset='test', categories=categories, shuffle=True, random_state=0, remove=remove)

## 完成
print (u"完成数据加载过程.耗时:%.3fs" % (time() - t_start))


# 获取加载的数据类别
data_train.target_names
data_train.target


### 获取加载数据的相关信息
def size_mb(docs):
    return sum(len(s.encode('utf-8')) for s in docs) / 1e6

categories = data_train.target_names
data_train_size_mb = size_mb(data_train.data)
data_test_size_mb = size_mb(data_test.data)

print (u'数据类型：', type(data_train.data))
print("%d文本数量 - %0.3fMB (训练数据集)" % (len(data_train.data), data_train_size_mb))
print("%d文本数量 - %0.3fMB (测试数据集)" % (len(data_test.data), data_test_size_mb))
print (u'训练集和测试集使用的%d个类别的名称：' % len(categories))
print(categories)


### 数据重命名
x_train = data_train.data
y_train = data_train.target
x_test = data_test.data
y_test = data_test.target


### 输出前5个样本
print (u' -- 前5个文本 -- ')
for i in range(5):
    print (u'文本%d(属于类别 - %s)：' % (i+1, categories[y_train[i]]))
    print (x_train[i])
    print ('\n\n')


### 文档转换为向量
## 转换
# 要求输入的数据类型必须是list，list中的每条数据是一个文本对象，单词是以空格分割开的
vectorizer = TfidfVectorizer(input='content', stop_words='english', max_df=0.5, sublinear_tf=True)
# 从训练集文档中找出作为特征属性的单词
x_train = vectorizer.fit_transform(data_train.data)  # x_train是稀疏的，scipy.sparse.csr.csr_matrix
x_test = vectorizer.transform(data_test.data)
print (u'训练集样本个数：%d，特征个数：%d' % x_train.shape)
print (u'停止词:\n')
print(vectorizer.get_stop_words())
## 获取最终的特征属性名称
print("最终的特征名称")
feature_names = np.asarray(vectorizer.get_feature_names())
print(feature_names.shape)
print(feature_names[:10])
print("转换之后的数据形式：")
print(x_train[0])


## 特征选择
ch2 = SelectKBest(chi2, k=1000)
x_train = ch2.fit_transform(x_train, y_train)
x_test = ch2.transform(x_test)
feature_names = [feature_names[i] for i in ch2.get_support(indices=True)]


### 使用不同的分类器对数据进行比较
print (u'分类器的比较：\n')
clfs = [
    [RidgeClassifier(), 'Ridge'],
    [KNeighborsClassifier(), 'KNN'],
    [MultinomialNB(), 'MultinomialNB'],
    [BernoulliNB(), 'BernoulliNB'],
    [RandomForestClassifier(n_estimators=200), 'RandomForest'],
    [SVC(), 'SVM'],
    [LinearSVC(loss='squared_hinge', penalty='l1', dual=False, tol=1e-4), 'LinearSVC-l1'],
    [LinearSVC(loss='squared_hinge', penalty='l2', dual=False, tol=1e-4), 'LinearSVC-l2']
]

## 将训练数据保存到一个列表中
result = []
for clf,name in clfs:
    # 计算算法结果
    a = benchmark(clf,name)
    # 追加到一个列表中，方便进行展示操作
    result.append(a)
    print ('\n')
## 将列表转换为数组
result = np.array(result)


### 获取需要画图的数据
result = [[x[i] for x in result] for i in range(5)]
training_time, test_time, training_err, test_err, clf_names = result

training_time = np.array(training_time).astype(np.float)
test_time = np.array(test_time).astype(np.float)
training_err = np.array(training_err).astype(np.float)
test_err = np.array(test_err).astype(np.float)


### 画图
x = np.arange(len(training_time))
plt.figure(figsize=(10, 7), facecolor='w')
ax = plt.axes()
b0 = ax.bar(x+0.1, training_err, width=0.2, color='#77E0A0')
b1 = ax.bar(x+0.3, test_err, width=0.2, color='#8800FF')
ax2 = ax.twinx()
b2 = ax2.bar(x+0.5, training_time, width=0.2, color='#FFA0A0')
b3 = ax2.bar(x+0.7, test_time, width=0.2, color='#FF8080')
plt.xticks(x+0.5, clf_names)
plt.legend([b0[0], b1[0], b2[0], b3[0]], (u'训练集错误率', u'测试集错误率', u'训练时间', u'测试时间'), loc='upper left', shadow=True)
plt.title(u'新闻组文本数据分类及不同分类器效果比较', fontsize=18)
plt.xlabel(u'分类器名称')
plt.grid(True)
plt.tight_layout(2)
plt.show()


algo = GaussianNB()
algo.fit(x_train.toarray(), y_train)
print("训练集效果:{}".format(algo.score(x_train.toarray(), y_train)))
print("测试集效果:{}".format(algo.score(x_test.toarray(), y_test)))

贝叶斯网络

把某个研究系统中涉及到的随机变量，根据是否条件独立绘制在一个有向图中，就形成了贝叶斯网络。

贝叶斯网络(Bayesian Network)，又称有向无环图模型(directed acyclic graphical model, DAG)，是一种概率图模型，根据概率图的拓扑结构，考察一组随机变量{X1 ,X2 ,...,Xn}及其N组条件概率分布(Conditional Probabililty Distributions, CPD)的性质。

当多个特征属性之间存在着某种相关关系的时候，使用朴素贝叶斯算法就没法解决这类问题，那么贝叶斯网络就是解决这类应用场景的一个非常好的算法。

一般而言，贝叶斯网络的有向无环图中的节点表示随机变量，可以是可观察到的变量，或隐变量，未知参数等等。连接两个节点之间的箭头代表两个随机变量之间的因果关系(也就是这两个随机变量之间非条件独立)，如果两个节点间以一个单箭头连接在一起，表示其中一个节点是“因” ，另外一个是“果” ，从而两节点之间就会产生一个条件概率值。

贝叶斯网络的关键方法是图模型，构建一个图模型我们需要把具有因果联系的各 个变量用箭头连在一起。贝叶斯网络的有向无环图中的节点表示随机变量。连接 两个节点的箭头代表此两个随机变量是具有因果关系的。 贝叶斯网络是模拟人的认知思维推理模式的，用一组条件概率以及有向无环图对 不确定性因果推理关系建模。

最简单的一个贝叶斯网络

全连接贝叶斯网络

每一对节点之间都有边连接：

“正常”贝叶斯网络

实际贝叶斯网络：判断是否下雨

有一天早晨，Bruce离开他的房子的时候发现他家花园中的草地是湿的，有两种可能，第一：昨天晚上下雨了，第二：他昨天晚上忘记关掉花园中的喷水器，接下来，他观察他的邻居Joe，发现他家花园中的草地也是湿的，因此，他推断，他家的草地湿了是因为昨天晚上下雨的缘故。

贝叶斯网络判定条件独立-01

在C给定的条件下，a和b被阻断(blocked)是独立的 条件独立：tail - to -tail

在C给定的条件下，a和b被阻断(blocked)是独立的。

• 条件独立：head- to -tail

贝叶斯网络判定条件独立-03