目录
3、逻辑回归(Lojit)与Fisher线性判别分析(LDA)如何选择?
(2)使用corr()方法计算特征变量间的相关系数矩阵,并绘制相关系数热力图
(5)将iris数据集随机分为训练集(70%)和测试集(30%),进行判别分析
1、背景知识
判别分析的目的就是从现有已知类别的样本数据中训练出一个判别函数,以后再有未知类别的样本进入,就利用建立的判别函数结合判别规则来判别其类别。各类判别问题的前提有所不同,进行划分或寻找判别函数的准则也可以不同,判别分析的方法有:距离判别法、Fisher判别法、贝叶斯判别法、逐步判别法,其中最常用的判别分析方法当属Fisher判别法。更详细的内容,可以查阅何晓群编著的《多元统计分析》。
2、Fisher判别法
- Fisher判别是根据方差分析的思想建立起来的,即按类内方差尽量小,类间方差尽量大的准则来求判别函数的
- Fisher判别的基本思想是投影,即按原来P维空间的自变量组合投影到维度较低的D维空间中去,然后在D维空间中再进行分类。投影的原则是使得每一类的离差尽可能小,而不同类间投影的离差尽可能大
- 对总体的分别不做任何的要求
3、逻辑回归(Lojit)与Fisher线性判别分析(LDA)如何选择?
虽然逻辑回归(Lojit)与Fisher线性判别分析(LDA)二者都可以解决分类问题,但一般来说,对多个类别进行预测时,LDA模型的效果要远高于逻辑回归模型。若一组预测变量需要被分成两类时,则一般使用逻辑回归。
4、使用线性判别分析建模前的数据预处理
响应变量必须是分类变量,特征变量则为数值型变量。
5、Fisher线性判别分析的应用领域举例
(1)在经济学中的应用:中小企业的破产模型
为了研究中小企业的破产模型,选定4个经济指标:X1为总负债率(现金收益/总负债)、X2为收益性指标(纯收入/总财产)、X3为短期支付能力(纯流动资产/流动负债)、X4为生产效率性指标(流动资产/纯销售额)。
- 对17个破产企业(第1类型)和21个正常运行企业(第2类型)进行了调查,得到关于上述4个指标的资料
- 现有8个未知类型的企业的4个经济指标数据,判断其属于破产企业一类还是属于正常企业一类
(2)在天气预报中应用
假如根据经验,下面两个因素是预报明天下雨还是不下雨的决定因素:
- 今天与昨天的湿度差,记为X1
- 今天的压温差(气压与温度之差),记为X2
今测得X1和X2,应预测明天下雨还是不下雨?
(3)医学领域
在医学实验中,研究者会经常使用LDA模型,去预测一组异常细胞,来判别它是否会对人体造成影响。
(4)生态领域
生态研究者会利用LDA模型预测环境的受污染程度:非常严重、严重、正常等,其特征变量包括硫 化物含量、温度、湿度等。
6、Fisher线性判别分析案例
本案例采用分类问题中的经典数据集:iris数据集,逐步介绍线性判别分析的python操作。数据可在https://faculty.sdu.edu.cn/qiang2chen2/zh_CN/jxzy/545776/content/1793.htm#jxzy中下载。
(1)导入相关模块,并查看数据
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
from sklearn.discriminant_analysis import QuadraticDiscriminantAnalysis
from sklearn.metrics import confusion_matrix
from sklearn.metrics import classification_report
from sklearn.metrics import cohen_kappa_score
iris = load_iris() #加载数据dir(iris)#查看iris数据集组成部分
输出结果:
['DESCR', 'data', 'data_module', 'feature_names', 'filename', 'frame', 'target', 'target_names'] #DESCR为数据集描述 #data为数据 #data_module为数据模块 #feature_names为特征变量名称 #filename为文件名 #frame为数据结构 #target为响应变量 #target_names为响应变量名称 iris.feature_names #查看特征变量x的名称输出结果:
['sepal length (cm)', 'sepal width (cm)', 'petal length (cm)', 'petal width (cm)']
结果显示,特征变量名称较为繁琐,带有空格和括号,故将特征变量的名称进行更新:
iris.feature_names = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width'] #重命名
X = pd.DataFrame(iris.data, columns=iris.feature_names) #变为数据框
X输出结果:
sepal_length sepal_width petal_length petal_width 0 5.1 3.5 1.4 0.2 1 4.9 3.0 1.4 0.2 2 4.7 3.2 1.3 0.2 3 4.6 3.1 1.5 0.2 4 5.0 3.6 1.4 0.2 ... ... ... ... ... 145 6.7 3.0 5.2 2.3 146 6.3 2.5 5.0 1.9 147 6.5 3.0 5.2 2.0 148 6.2 3.4 5.4 2.3 149 5.9 3.0 5.1 1.8 150 rows × 4 columns
(2)使用corr()方法计算特征变量间的相关系数矩阵,并绘制相关系数热力图
使用seaborn模块的heatmap()函数将相关系数矩阵进行可视化,通过
X.corr()#计算pearson相关系数
sns.heatmap(X.corr(), cmap='RdBu_r', annot=True)#annot=True表示显示数值注释
plt.savefig("squares3.png",
bbox_inches ="tight",
pad_inches = 1,
transparent = True,
facecolor ="w",
edgecolor ='w',
dpi=300,
orientation ='landscape')#保存高清图片输出结果:
(3)查看响应变量取值
y = iris.target
y输出结果:
array([0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2])
(4)使用全样本(不区分训练集和测试集)进行线性判别分析
dir(LinearDiscriminantAnalysis())#查看用类的属性
输出结果:
['__annotations__', '__class__', '__delattr__', '__dict__', '__dir__', '__doc__', '__eq__', '__format__', '__ge__', '__getattribute__', '__getstate__', '__gt__', '__hash__', '__init__', '__init_subclass__', '__le__', '__lt__', '__module__', '__ne__', '__new__', '__reduce__', '__reduce_ex__', '__repr__', '__setattr__', '__setstate__', '__sizeof__', '__str__', '__subclasshook__', '__weakref__', '_check_feature_names', '_check_n_features', '_estimator_type', '_get_param_names', '_get_tags', '_more_tags', '_parameter_constraints', '_predict_proba_lr', '_repr_html_', '_repr_html_inner', '_repr_mimebundle_', '_sklearn_auto_wrap_output_keys', '_solve_eigen', '_solve_lstsq', '_solve_svd', '_validate_data', '_validate_params', 'covariance_estimator', 'decision_function', 'fit', 'fit_transform', 'get_feature_names_out', 'get_params', 'n_components', 'predict', 'predict_log_proba', 'predict_proba', 'priors', 'score', 'set_output', 'set_params', 'shrinkage', 'solver', 'store_covariance', 'tol', 'transform']
model = LinearDiscriminantAnalysis()#创建线性判别分析的一个实例
model.fit(X, y) #模型拟合
model.score(X, y) #计算预测准确率
输出结果:0.98model.priors_#查看各类别的先验概率
输出结果:array([0.33333333, 0.33333333, 0.33333333])model.means_#查看4个特征变量的分组均值
输出结果:
array([[5.006, 3.428, 1.462, 0.246], [5.936, 2.77 , 4.26 , 1.326], [6.588, 2.974, 5.552, 2.026]])
model.explained_variance_ratio_#查看线性判元对于组间方差的贡献输出结果:
array([0.9912126, 0.0087874])model.scalings_#查看线性判元的估计系数
lda_loadings = pd.DataFrame(model.scalings_, index=iris.feature_names, columns=['LD1', 'LD2'])#将线性判别元系数变为数据框展示
lda_loadings输出结果:
LD1 LD2 sepal_length 0.829378 0.024102 sepal_width 1.534473 2.164521 petal_length -2.201212 -0.931921 petal_width -2.810460 2.839188
其中第1列为第1个线性判别元LD1对于4个特征变量的线性判别载荷,而第2列为第2个线性判别元LD2对于4个特征变量的线性判别载荷。针对已估计的线性判别模型,使用transform()方法,得到线性判别得分。
lda_scores = model.fit(X, y).transform(X) #或lda_scores = model.fit_transform(X, y)
lda_scores.shape#查看线性判别得分数据形状输出结果:
(150, 2)lda_scores[:5, :]#展示前5行观测值
输出结果:
array([[ 8.06179978, 0.30042062], [ 7.12868772, -0.78666043], [ 7.48982797, -0.26538449], [ 6.81320057, -0.67063107], [ 8.13230933, 0.51446253]])LDA_scores = pd.DataFrame(lda_scores, columns=['LD1', 'LD2'])#将lda_scores数据由数组变为数据框,并对每一列命名
LDA_scores['Species'] = iris.target#在数据框中加上响应变量
LDA_scores.head()输出结果:
LD1 LD2 Species 0 8.061800 0.300421 0 1 7.128688 -0.786660 0 2 7.489828 -0.265384 0 3 6.813201 -0.670631 0 4 8.132309 0.514463 0 ... ... ... ... 145 -5.645003 1.677717 2 146 -5.179565 -0.363475 2 147 -4.967741 0.821141 2 148 -5.886145 2.345091 2 149 -4.683154 0.332034 2 150 rows × 3 columns
d = {0: '山鸢尾花', 1: '杂色鸢尾花', 2: '维吉尼亚鸢尾花'}
LDA_scores['Species'] = LDA_scores['Species'].map(d) #将变量Species的取值从{0,1,2}变为{山鸢尾花,杂色鸢尾花,维吉尼亚鸢尾花}
LDA_scores输出结果
LD1 LD2 Species 0 8.061800 0.300421 山鸢尾花 1 7.128688 -0.786660 山鸢尾花 2 7.489828 -0.265384 山鸢尾花 3 6.813201 -0.670631 山鸢尾花 4 8.132309 0.514463 山鸢尾花 ... ... ... ... 145 -5.645003 1.677717 维吉尼亚鸢尾花 146 -5.179565 -0.363475 维吉尼亚鸢尾花 147 -4.967741 0.821141 维吉尼亚鸢尾花 148 -5.886145 2.345091 维吉尼亚鸢尾花 149 -4.683154 0.332034 维吉尼亚鸢尾花 150 rows × 3 columns
import matplotlib
import matplotlib.pyplot as plt
matplotlib.rcParams['font.sans-serif'] = ['FangSong'] # 用仿宋字体显示中文
plt.rcParams['axes.unicode_minus'] = False # 显示负数不乱码
sns.scatterplot(x='LD1', y='LD2', data=LDA_scores, hue='Species')
plt.savefig("squares4.png",
bbox_inches ="tight",
pad_inches = 1,
transparent = True,
facecolor ="w",
edgecolor ='w',
dpi=300,
orientation ='landscape')输出结果:
可以看出LD1这个维度上很好区分不同的鸢尾花品种。
观察每个特征的分布情况:
from sklearn import datasets
import matplotlib.pyplot as plt
import numpy as np
import math
# prepare the data
iris = datasets.load_iris()X = iris.data
y = iris.target
iris.feature_names = ['sepal_length', 'sepal_width', 'petal_length', 'petal_width']
names = iris.feature_names
iris.target_names= ['山鸢尾花', '杂色鸢尾花', '维吉尼亚鸢尾花']
labels = iris.target_names
y_c = np.unique(y)
fig, axes = plt.subplots(2, 2, figsize=(12, 6))
for ax, column in zip(axes.ravel(), range(X.shape[1])):
# set bin sizes
min_b = math.floor(np.min(X[:, column]))
max_b = math.ceil(np.max(X[:, column]))
bins = np.linspace(min_b, max_b, 25)
# plotting the histograms
for i, color in zip(y_c, ('gold', 'cyan', 'violet')):
ax.hist(X[y == i, column], color=color, label='class %s' % labels[i],
bins=bins, alpha=0.5, )
ylims = ax.get_ylim()
# plot annotation
l = ax.legend(loc='upper right', fancybox=True, fontsize=8)
l.get_frame().set_alpha(0.5)
ax.set_ylim([0, max(ylims) + 2])
ax.set_xlabel(names[column])
ax.set_title('Iris histogram feature %s' % str(column + 1))
# hide axis ticks
ax.tick_params(axis='both', which='both', bottom=False, top=False, left=False, right=False,
labelbottom=True, labelleft=True)
# 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)
axes[0][0].set_ylabel('count')
axes[1][0].set_ylabel('count')
fig.tight_layout()plt.savefig("squares5.png",
bbox_inches ="tight",
pad_inches = 1,
transparent = True,
facecolor ="w",
edgecolor ='w',
dpi=300,
orientation ='landscape')输出结果:
通过简单的特征图形表示,petal_length和petal_width更适合分类。下面只采用这两个特征变量进行线性判别分析,直观的画出线性判别的决策边界:
import numpy as np
X = pd.DataFrame(iris.data, columns=iris.feature_names) #变为数据框
X2 = X.iloc[:, 2:4]
model = LinearDiscriminantAnalysis()
model.fit(X2, y)
model.score(X2, y)输出结果:
0.96model.explained_variance_ratio_#考察线性判别元对组间方差的贡献率
输出结果;
array([0.99470499, 0.00529501])
结果显示,即使只用两个 特征变量,预测准确率仍可达到96%,而第1个线性判别元即可解释约99.47%的组间方差。
from mlxtend.plotting import plot_decision_regions#导入模块
plot_decision_regions(np.array(X2), y, model)#画出线性判别决策边界,plot_decision_regions要求数据矩阵为数组
plt.xlabel('petal_length')
plt.ylabel('petal_width')
plt.title('Decision Boundary for LDA')
plt.savefig("squares6.png",
bbox_inches ="tight",
pad_inches = 1,
transparent = True,
facecolor ="w",
edgecolor ='w',
dpi=300,
orientation ='landscape')输出结果:
通过图像可以看出,线性判别分析的决策边界确实为线性,而且可以较好地区分三种鸢尾花。
(5)将iris数据集随机分为训练集(70%)和测试集(30%),进行判别分析
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, stratify=y, random_state=123)
model = QuadraticDiscriminantAnalysis()
model.fit(X_train, y_train)
model.score(X_test, y_test)#测试集的预测准确率输出结果:
0.9777777777777777prob = model.predict_proba(X_test)#计算预测测试集中每个样品的分类概率
prob[:3]
输出结果:array([[1.30728298e-117, 4.25364516e-001, 5.74635484e-001], [4.80474602e-185, 1.17558005e-003, 9.98824420e-001], [1.12275433e-056, 9.99994095e-001, 5.90465247e-006]])pred = model.predict(X_test)#实现训练的模型预测测试集中每个样品的种类
pred[:5]输出结果:
array([2, 2, 1, 1, 2])confusion_matrix(y_test, pred)#计算混淆矩阵
输出结果:
array([[15, 0, 0], [ 0, 14, 1], [ 0, 0, 15]], dtype=int64)print(classification_report(y_test, pred)) #打印基于混淆矩阵的一系列预测效果指标
cohen_kappa_score(y_test, pred)#计算科恩kappa指标输出结果:
precision recall f1-score support 0 1.00 1.00 1.00 15 1 1.00 0.93 0.97 15 2 0.94 1.00 0.97 15 accuracy 0.98 45 macro avg 0.98 0.98 0.98 45 weighted avg 0.98 0.98 0.98 45 0.9666666666666667
结果显示,科恩kappa值约为0.967,说明预测值和真实值之间非常一致,线性判别分析对于数据集iris达到了很好的分类效果。
(6)二次判别分析
同样选取两个特征变量进画出其决策边界:
X2 = X.iloc[:, 2:4]
model = QuadraticDiscriminantAnalysis()#二次判别分析
model.fit(X2, y)
model.score(X2, y)
plot_decision_regions(np.array(X2), y, model)
plt.xlabel('petal_length')
plt.ylabel('petal_width')
plt.title('Decision Boundary for QDA')
plt.savefig("squares7.png",
bbox_inches ="tight",
pad_inches = 1,
transparent = True,
facecolor ="w",
edgecolor ='w',
dpi=300,
orientation ='landscape')输出结果:
结果显示,二次判别分析的决策边界为二次函数(抛物线),也能较好地区分三种鸢尾花。
更多优质内容持续发布中,请移步主页查看。
点赞+关注,下次不迷路!