LDA（线性判别式分析）实践与学习

LDA实验报告

• 实验概述

LDA，线性判别式分析（Linear Discriminant Analysis），是模式识别的经典算法。基本思想是将高维的模式样本投影到最佳鉴别矢量空间，以达到抽取分类信息和压缩特征空间维数的效果，投影后类内方差最小，类间方差最大。

• 实验与思考

1. 关于LDA降维效果的思考

左图是单纯最大化类间距离可是并不能把红蓝两类很好的区分开来，注意这时候w和(m1−m2)是平行的。

右图则是使用LDA降维的结果。

四维降到一维

四维降到二维

 

1. LDA直线的斜率为最大的特征值对应的特征向量w的转置，而截距随意，垂直于LDA的直线’ 能够将两堆点进行分类，截距通过各堆点与垂足的直线截距计算出范围，然后便可以随机生成一个适合的截距。

 

1. LDA不适合对非正态分布样本进行降维。

     

      上图中红色区域表示一类样本，蓝色区域表示另一类，由于是2类，所以最多投影到1维上。不管在直线上怎么投影，都难使红色点和蓝色点内部凝聚，类间分离。如果样本数据由独立非高斯分布的隐含因子产生，隐含因子个数等于特征数，比如说信号这样的数据，就可以用ICA。

      实验代码

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import random
# 计算均值,要求输入数据为numpy的矩阵格式，行表示样本数，列表示特征


def meanX(data):
    return np.mean(data, axis=0)
    # axis=0表示按照列来求均值，如果输入list,则axis=1

# 计算类内离散度矩阵子项si
def compute_si(xi):
    n = xi.shape[0]
    ui = meanX(xi)
    si = 0
    for i in range(0, n):
        si = si + (xi[i, :] - ui).T * (xi[i, :] - ui)
    return si

# 计算类间离散度矩阵Sb
def compute_Sb(x1, x2):
    dataX = np.vstack((x1, x2))
    u1 = meanX(x1)
    u2 = meanX(x2)
    u = meanX(dataX)
    Sb = (u-u1).T * (u-u1) + (u-u2).T * (u-u2)
    return Sb


def LDA(x1, x2):
    s1 = compute_si(x1)
    s2 = compute_si(x2)
    n1 = len(s1)
    n2 = len(s2)
    Sw = (n1*s1+n2*s2)/(n1+n2)
    Sb = compute_Sb(x1, x2)

    # 求最大特征值对应的特征向量
    eig_value, vec = np.linalg.eig(np.mat(Sw).I * Sb)  # 特征值和特征向量
    index_vec = np.argsort(-eig_value)  # 对eig_value从大到小排序，返回索引
    eig_index = index_vec[:1]  # 取出最大的特征值的索引
    w = vec[:, eig_index]  # 取出最大的特征值对应的特征向量
    return w


def createDataSet():
    X1 = np.mat(np.random.random((10, 2)) * 10 + 15)  # 类别A
    X2 = np.mat(np.random.random((10, 2)) * 10 + 2)  # 类别B
    return X1, X2


def plotFig(data1, data2, k):
    fig = plt.figure(figsize=(5, 5))

    plt.ylim(0, 30)
    plt.xlim(0, 30)
    ax = fig.add_subplot(111)
    x1 = data1[0].tolist()[0]
    y1 = data1[1].tolist()[0]
    x2 = data2[0].tolist()[0]
    y2 = data2[1].tolist()[0]
    ax.scatter(x=x1, y=y1, color='red')
    ax.scatter(x=x2, y=y2, color='blue')
    # 画LDA直线
    x = np.linspace(0, 30, 300)
    line1 = k*x
    plt.plot(x, line1, color="yellow", linewidth=2)

    # 画投影点
    m = k
    b = 0
    minb = 0
    maxb = 100
    xy = data1.T.tolist()
    for x0, y0 in xy:
        xx = (m*y0+x0-m*b)/(m**2+1)
        yy = (m**2*y0+m*x0+b)/(m**2+1)
        ax.scatter(xx, yy, color='r')
        plt.plot([x0, xx], [y0, yy], color='m')
        maxb = min(maxb, y0-(y0-yy)/(x0-xx)*x0)

    xy = data2.T.tolist()
    for x0, y0 in xy:
        xx = (m*y0+x0-m*b)/(m**2+1)
        yy = (m**2*y0+m*x0+b)/(m**2+1)
        ax.scatter(xx, yy,  color='b')
        plt.plot([x0, xx], [y0, yy],  color='c')
        minb = max(minb, y0-(y0-yy)/(x0-xx)*x0)

    # 画与LDA直线垂直的直线
    line2 = (-1/k)*x+random.uniform(minb, maxb)
    plt.plot(x, line2, color="green", linewidth=2)

    plt.show()


if __name__ == "__main__":
    x1, x2 = createDataSet()
    w = LDA(x1, x2)
    print('w', w)
    plotFig(x1.T, x2.T, w[0, 0]/w[1, 0])

# -*- coding: utf-8 -*-
import numpy as np
from matplotlib import pyplot as plt
from sklearn import datasets
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA


def read_iris():
    from sklearn.datasets import load_iris
    from sklearn import preprocessing
    data_set = load_iris()
    data_x = data_set.data
    label = data_set.target + 1
    preprocessing.scale(data_x, axis=0, with_mean=True,
                        with_std=True, copy=False)
    return data_x, label
 # 特征均值,计算每类的均值，返回一个向量


def class_mean(data, label, clusters):
    mean_vectors = []
    for cl in range(1, clusters+1):
        mean_vectors.append(np.mean(data[label == cl, ], axis=0))
    # print mean_vectors
    return mean_vectors

# 计算类内散度


def within_class_SW(data, label, clusters):
    m = data.shape[1]
    S_W = np.zeros((m, m))
    mean_vectors = class_mean(data, label, clusters)
    for cl, mv in zip(range(1, clusters+1), mean_vectors):
        class_sc_mat = np.zeros((m, m))
        # 对每个样本数据进行矩阵乘法
        for row in data[label == cl]:
            row, mv = row.reshape(4, 1), mv.reshape(4, 1)
            class_sc_mat += (row-mv).dot((row-mv).T)
        S_W += class_sc_mat
    return S_W


def between_class_SB(data, label, clusters):
    m = data.shape[1]
    all_mean = np.mean(data, axis=0)
    S_B = np.zeros((m, m))
    mean_vectors = class_mean(data, label, clusters)
    for cl, mean_vec in enumerate(mean_vectors):
        n = data[label == cl+1, :].shape[0]
        mean_vec = mean_vec.reshape(4, 1)  # make column vector
        all_mean = all_mean.reshape(4, 1)  # make column vector
        S_B += n * (mean_vec - all_mean).dot((mean_vec - all_mean).T)
    return S_B


def lda():
    data, label = read_iris()
    clusters = 3
    S_W = within_class_SW(data, label, clusters)
    S_B = between_class_SB(data, label, clusters)
    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)))
    print('Matrix W:\n', W.real)
    print(data.dot(W))
    return W


def D2_plot_lda():
    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    target_names = iris.target_names
    lda = LDA(n_components=2)
    Y = lda.fit(X, y).transform(X)
    data, labels = read_iris()
    ax = plt.subplot(111)
    for label, marker, color in zip(range(1, 4), ('^', 's', 'o'), ('red', 'yellow', 'blue')):
        plt.scatter(x=Y[:, 0][labels == label],
                    y=Y[:, 1][labels == label],
                    marker=marker,
                    color=color,
                    alpha=0.5,
                    )
    plt.grid()
    plt.show()


def D1_plot_lda():
    iris = datasets.load_iris()
    X = iris.data
    y = iris.target
    lda = LDA(n_components=1)
    X_r2 = lda.fit(X, y).transform(X)
    X_Zreo = np.zeros(X_r2.shape)
    for c, i in zip('ryb', [0, 1, 2]):
        plt.scatter(X_r2[y == i], X_Zreo[y == i], c=c)
    plt.grid()
    plt.show()


if __name__ == "__main__":
    lda()
    D2_plot_lda()
    D1_plot_lda()