『科学计算』高斯判别分析模型实现

和上一篇一样，本部分的理论建议自行学习cs229或者其他的高斯判别分析模型介绍文章。

1.模型简介

高斯判别分析模型是一种生成模型，而逻辑回归是一种判别模型，生成模型和判别模型的详细了解可参考这篇文章：

         http://blog.sciencenet.cn/home.php?mod=space&uid=248173&do=blog&id=227964

简单的来说，我们的目标都是p(y|x)，判别模型是构造一个函数f(x)去逼近p(y|x)，而对于生成模型则是通过贝叶斯公式p(y|x) = p(x|y)p(y)/p(x)，求得p(x|y)和p(y)来间接得到p(y|x)。

        

首先，高斯判别分析模型对变量x和y有如下假设：

          

这样，可以给出概率密度函数：

 

2.评价

         该模型的对数似然函数如下：

 

        

3.优化

         对各个参数进行求导后令等式为0，得到：

         

　　  Φ是训练样本中结果 y=1 占有的比例。
　　  μ0是 y=0 的样本中特征均值。
　　  μ1是 y=1 的样本中特征均值。
　　  Σ是样本特征方差均值。

代码如下，

import numpy as np
import pandas as pd
from sklearn.datasets import load_iris

# iris = pd.read_csv('http://aima.cs.berkeley.edu/data/iris.csv',
#                    names=['col0','col1','col2','col3','class'])
# dummy = pd.get_dummies(iris['col3'])
# iris = pd.concat([iris, dummy], axis=1)

iris = load_iris()
X = iris.data[:, 0:2]
Y = np.array(pd.get_dummies(iris.target)[0])
# Y = Y[Y[0]==1.]
# print(X[Y==0].mean(axis=0))

def GDA(X, Y):
    theta1 = Y.mean()
    theta0 = 1-theta1
    mu1 = X[Y==1].mean(axis=0)
    mu0 = X[Y==0].mean(axis=0)

    X1 = X[Y==1]
    X0 = X[Y==0]
    A = np.dot(X1.T, X1) - len(Y[Y==1])*np.dot(mu1.reshape(X.shape[1],1), mu1.reshape(X.shape[1],1).T)
    B = np.dot(X0.T, X0) - len(Y[Y==0])*np.dot(mu0.reshape(X.shape[1],1), mu0.reshape(X.shape[1],1).T)
    sigma = (A+B)/len(X)

    return theta1, mu1, mu0, sigma

if __name__=='__main__':
    theta1, mu1, mu0, sigma = GDA(X, Y)
    print(theta1,
          '
', mu1,
          '
', mu0,
          '
', sigma)

我们来检查一下数据，

X.shape
Out[2]:
(150, 2)

Y.shape
Out[3]:
(150,)

由于是二分类问题，实际上我们Y的one_hot只表示属于类别1（1）和其他类别（2）两种标签。

实际上iris是有4个特征的，我们只取了前两个，为什么呢。。。因为我想可视化，高维特征不能可视化233，

简单的把输出

0.333333333333
[ 5.006 3.418]
[ 6.262 2.872]
[[ 0.33055867 0.113388 ]
[ 0.113388 0.12050267]]

导入『科学计算』可视化二元正态分布一节中的可视化函数即可，

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import axes3d
from matplotlib import cm
import matplotlib as mpl

num = 500
l = np.linspace(0,10,num)
X, Y =np.meshgrid(l, l)
pos = np.concatenate((np.expand_dims(X,axis=2),np.expand_dims(Y,axis=2)),axis=2)

u1 = np.array([5.006, 3.418])
o1 = 3*np.array([[0.33055867, 0.113388],
                 [0.113388, 0.12050267]])
a1 = (pos-u1).dot(np.linalg.inv(o1))
b1 = np.expand_dims(pos-u1,axis=3)
Z1 = np.zeros((num,num), dtype=np.float32)

u2 = np.array([6.262, 2.872])
o2 = 3*np.array([[0.33055867, 0.113388],
                 [0.113388, 0.12050267]])
a2 = (pos-u2).dot(np.linalg.inv(o2))
b2 = np.expand_dims(pos-u2,axis=3)
Z2 = np.zeros((num,num), dtype=np.float32)

for i in range(num):
    Z1[i] = [np.dot(a1[i,j],b1[i,j]) for j in range(num)]
    Z2[i] = [np.dot(a2[i,j],b2[i,j]) for j in range(num)]
Z1 = np.exp(Z1*(-0.5))/(2*np.pi*np.linalg.det(o1))
Z2 = np.exp(Z2*(-0.5))/(2*np.pi*np.linalg.det(o1))

Z = Z1 + Z2

fig = plt.figure()
ax = fig.add_subplot(211,projection='3d')
ax.plot_surface(X, Y, Z, rstride=5, cstride=5, alpha=0.5, cmap=mpl.cm.rainbow)

ax.contour(X,Y,Z1,10,zdir='z',offset=0,cmap=cm.coolwarm)
ax.contour(X,Y,Z2,10,zdir='z',offset=0,cmap=cm.coolwarm)
ax.contour(X, Y, Z, zdir='x', offset=-0,cmap=mpl.cm.winter)
ax.contour(X, Y, Z, zdir='y', offset= 10,cmap= mpl.cm.winter)
'''
mpl.cm.rainbow
mpl.cm.winter
mpl.cm.bwr  # 蓝，白，红
cm.coolwarm
'''

ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('Z')
plt.show()

ax2 = fig.add_subplot(212)
cs = ax2.contour(X,Y,Z1)
ax2.clabel(cs, inline=1, fontsize=20)
cs2 = ax2.contour(X,Y,Z2)
ax2.clabel(cs2, inline=1, fontsize=20)

输出图像如下（调整了一下坐标显示，要不然显示不全），

换了个颜色233，