LDA线性判别分析

LAD概括

LAD主要用于在高维数据的分类中,将数据按照线性模型进行降维;使得数据得到很好的分类同时避免数据过度拟合,与PAC不同的是

LAD主要用于有监督的学习,用于对目标分类;约束条件则是类间的协方差最大,同类内的协方差最小化
PAC则是无监督的学习,将样本投影到几个正交的方向,同时同类内的方差越大越好,尽可能的保留样本的全部信息


两类的LAD问题可以看作是把所有样本都投影到一个方向上,然后在这个一维空间中确定一个分类的阈值。过这个阈值点且与投影方向垂直的超平面就是两类的分类面。

公式推导

  • 映射直线
    y=wTx
  • 约束条件则是类间的协方差最大,同类内的协方差最小化

J=||wTμ0−wTμ1||22/wT(∑0+∑1)w=wT(μ0−μ1)(μ0−μ1)Tw/wT(∑0+∑1)w

这里涉及到欧氏距离也就是矩阵的平方等于 原矩阵*矩阵的转置

  • 目的使得约束条件J最大化
    类间协方差矩阵 Sb=(μ0−μ1)(μ0−μ1)T
    类间协方差矩阵 Sw=∑x∈D0(x−μ0)(x−u0)T+∑x∈D1(x−μ1)(x−u1)T
    重写J: J=wTSbw/wTSww

  • 关于W的确定,由于分子分母是w的平方,已经将w实数化;因此方程的解与w的大小已经无关,仅仅是方向上有关了

拉格朗日乘子法

令wTSww=1
c(w)=wTSbw−λ(wTSww−1)
涉及到矩阵的求导
dc/dw=2Sbw−2λSww=0
Sbw=λSww
SbW 类间散度矩阵 方向是固定的U0-U1
w=Sw -1(u0-u1)


基于两个以上特征值的分类

原理是一样的,在计算类间散度矩阵是多了两个参数
SB=∑i=1-c Ni(mmi−mm)(mmi−mm)T
mm 是全局均值,而 mmi 和 Ni 是每类样本的均值和样本数。

求出特征值和特征序列

矩阵 S−1WSB =W
按照特征值排序,将原来多维度的空间映射到对应的空间

代码实现

数据获取与格式整理

#coding=UTF-8
import numpy as np
import pandas as pd 
from sklearn.preprocessing import LabelEncoder
from matplotlib import pyplot as plt
####构造索引与对应的名称关系
feature_dict = {i:label for i,label in zip(
                range(4),
                  ('sepal length in cm',
                  'sepal width in cm',
                  'petal length in cm',
                  'petal width in cm', ))}
'''
{0: 'sepal length in cm',
 1: 'sepal width in cm', 
 2: 'petal length in cm', 
 3: 'petal width in cm'}
'''
##读取数据
df=pd.io.parsers.read_csv(
    filepath_or_buffer='https://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data',
    header=None,
    sep=',',
    )
###定义列名
df.columns=[l for i,l in sorted(feature_dict.items())]+["classLab"]
#删除文件末尾的空行
df.dropna(how="all",inplace=True)
###看一下最后10行时啥样
#print(df.tail())
'''
     sepal length in cm  sepal width in cm  petal length in cm  \
145                 6.7                3.0                 5.2   
146                 6.3                2.5                 5.0   
147                 6.5                3.0                 5.2   
148                 6.2                3.4                 5.4   
149                 5.9                3.0                 5.1   

     petal width in cm        classLab  
145                2.3  Iris-virginica  
146                1.9  Iris-virginica  
147                2.0  Iris-virginica  
148                2.3  Iris-virginica  
149                1.8  Iris-virginica 
'''

将类别数字化

X=df[[0,1,2,3]].values #获取四种属性值
y=df['classLab'].values #获取类别名
enc = LabelEncoder()
label_encoder = enc.fit(y)
y = label_encoder.transform(y) + 1

label_dict = {1: 'Setosa', 2: 'Versicolor', 3:'Virginica'}
# print(len(y))
'''
数据处理完成

计算每类花对应的特征值的均值

##分别对三种花求其在4种属性上的均值
np.set_printoptions(precision=4)
mean_vectors=[]
for i in range(1,4):
    ##对4中属性150 X 4 的矩阵按列求均值
    mean_vectors.append(np.mean(X[y==i],axis=0))
    #print('mean Vector calss %s: %s\n' %(i,mean_vectors[i-1]))
'''
mean Vector calss 1: [5.006 3.418 1.464 0.244]

mean Vector calss 2: [5.936 2.77  4.26  1.326]

mean Vector calss 3: [6.588 2.974 5.552 2.026]
'''

类内散度矩阵

S_W=np.zeros((4,4))
for cl,mv in zip(range(1,4),mean_vectors):
    ####每个类的类内散度矩阵
    class_sc_mat=np.zeros((4,4))
    for row in X[y==cl]:
        ##把对应的值与特征值进行运算
        '''
        [1,第一个特征值
        2, 第2个特征值
        3, 第3个特征值
        4  第4个特征值
        ]

        '''
        row, mv = row.reshape(4,1), mv.reshape(4,1)
        class_sc_mat+=(row-mv).dot((row-mv).T)
    S_W+=class_sc_mat

类间散度矩阵

overall_mean=np.mean(X,axis=0)
S_b=np.zeros((4,4))
for i,men_vector in enumerate(mean_vectors):
    ###获取样本数目
    n=X[y==i+1,:].shape[0]
    men_vector=men_vector.reshape(4,1)
    overall_mean=overall_mean.reshape(4,1)
    S_b+=n*(men_vector-overall_mean).dot((men_vector- overall_mean).T)
print('between-class Scatter Matrix:\n', S_b)

获取特征值和特征向量

eig_vals, eig_vecs = np.linalg.eig(np.linalg.inv(S_W).dot(S_b))

for i in range(len(eig_vals)):
    eigvec_sc = eig_vecs[:,i].reshape(4,1)   
    print('\nEigenvector {}: \n{}'.format(i+1, eigvec_sc.real))
    print('Eigenvalue {:}: {:.2e}'.format(i+1, eig_vals[i].real))

根据特征值大小选取特征向量降维

eig_pairs = [(np.abs(eig_vals[i]), eig_vecs[:,i]) for i in range(len(eig_vals))]
eig_pairs = sorted(eig_pairs, key=lambda k: k[0], reverse=True)
W = np.hstack((eig_pairs[0][1].reshape(4,1), eig_pairs[1][1].reshape(4,1)))

映射到新的空间

X_lda = X.dot(W)
assert X_lda.shape == (150,2)

作图

def plot_step_lda():

    ax = plt.subplot(111)
    for label,marker,color in zip(
        range(1,4),('^', 's', 'o'),('blue', 'red', 'green')):

        plt.scatter(x=X_lda[:,0].real[y == label],
                y=X_lda[:,1].real[y == label],
                marker=marker,
                color=color,
                alpha=0.5,
                label=label_dict[label]
                )

    plt.xlabel('LD1')
    plt.ylabel('LD2')

    leg = plt.legend(loc='upper right', fancybox=True)
    leg.get_frame().set_alpha(0.5)
    plt.title('LDA: Iris projection onto the first 2 linear discriminants')

    # hide axis ticks
    plt.tick_params(axis="both", which="both", bottom="off", top="off",  
            labelbottom="on", left="off", right="off", labelleft="on")

    # remove axis spines
    ax.spines["top"].set_visible(False)  
    ax.spines["right"].set_visible(False)
    ax.spines["bottom"].set_visible(False)
    ax.spines["left"].set_visible(False)    

    plt.grid()
    plt.tight_layout
    plt.savefig("LDA.png")

plot_step_lda()

Reference

推荐看老外的英文介绍
老外的翻译版
机器学习-线性判别分析.周志华x
LAC
拉格朗日乘子法
PAC
二分类LAD