Loading...
墨滴

精益修身

2021/04/25  阅读:35  主题:嫩青

多元统计之判别分析

MV之判别分析

精益修身

1. 判别分析的理论

1.1 判别分析介绍

银行在判断申请者能否偿还贷款的时候,需要根据申请人的信息,将其进行分类。这个过程包括两个步骤:

第一:在银行既有的数据中,有能够偿还贷款的,也有不能偿还贷款的。第一步,我们要根据历史数据找出区分规则来区分这两类人。这一步叫做判别分析,有时也被称为描述性判别分析(descriptive discriminant analysis

第二:将新申请贷款的人的特征输入到判别规则中,对其进行分类。这一步叫做分类分析,有时也被称为预测性判别分析(predictive discriminant analysis

由于判别分析的数据包含有类别,所以判别分析属于有监督的学习(supervied learning,相反,聚类分析的数据没有类别,需要自己找到类别,属于无监督的学习(unsupervised learning

1.2 判别分析的区分方法

在Stata中,判别分析有四个不同的方法,分别是:

1、knn,k近邻判别分析;

2、lda,线性判别分析;

3、logistic判别分析;

4、qda,二次判别分析。

上述四种方法的差别在于确定组i的密度函数 是不同的。lda假定各组都是多元正态且协方差矩阵相等;qda同样假定多元正态,但允许协方差矩阵不等;logistic使用多元logistic模型来得到后验概率,knn使用非参数方法估计 。下面仅介绍knn和lda两种区分办法。

1.2.1 k近邻判别分析

knn是一种非参数判别方法,基于k个最近邻居的观测值。

1.2.1.1 原理
  1. 分组概率

如果数据有 组,第 组有 个观测值,第 组的先验概率为 表示某个观测,且有 区分变量(discriminating variables。虽然在Stata中,观测值都是横向的,但这里我们要把 看作是一个列向量 表示组 的密度函数, 表示在属于 的条件下观测到 的概率。 表示看到观测值 的条件下组 的后验概率。根据贝叶斯定理,我们有:

如果用 代替 ,则上式可以改写为:

对于knn判别方法,我们令 表示 个邻近的邻居中属于组 的个数。后验概率公式可改为:

  1. K的取值

    没有明确的定法,可以通过不同的取值,选择错分率(error rate)最低的K。有人建议是 ,即某一组样本数的平方根。

1.2.1.2 操作命令

语法结构:

discrim knn varlist [if] [in] [weight], group(groupvar) k(#) [options]

选项:

group(groupvar) 必选项,分组依据,括号内是分组变量。

k(#),必选项,在某一距离内,选择#个最近的观测。

priors(priors),设定先验概率,允许priors(equal)(默认),priors(proportional),priors(matname),pr iors(matrix exp)。

ties(ties),出现距离相同的情况叫tie(平局)。允许ties(missing)(默认),ties(random),ties(first),tie s(nearest)。

measure(measure),设定相似或相异的度量方法,默认是L2,即欧氏距离。

s2d(standard | oneminus),设定相似如何转化为相异,默认是standard。

mahalanobis,表示在计算相异前进行马氏变换。

1.2.1.3 估计后命令
estat classtable, 分类表
estat errorrate, 分错率估计
estat grsummarize, 组摘要
estat list, 分类列表
estat summarize, 估计样本摘要
estimates, 估计结果目录
predict, 预测分组和后验概率

1.2.2 Fisher线性判别法(LDA)

discrim lda 和 candisc(典型线性判别分析)的计算结果完全相同。

Fisher线性判别法假定数据服从正态分布,且协方差矩阵相同。Fisher(1936)通过寻找区分变量的线性组合来进行线性判别分析,该线性组合提供了组之间的最大分隔(最初是两个组,但后来扩展到多个组)。组的最大分隔由 的特征值和特征向量确定,其中 是组间平方和和叉积(SSCP)矩阵, 是组内SSCP矩阵。 的特征值和特征向量组成了Fisher线性判别函数。

线性判别函数是由特征向量构成,第一个线性判别函数所关联的特征值最大。该第一判别函数将原始变量线性转换到一个维度,在这个维度上,组均值之间有最大的分隔。第二大特征值关联的特征向量是第二线性判别函数。它提供了一个新的维度,这一维度与第一维度不相关(但通常不正交)。第二个判别函数在第二维中提供了最大的组间隔。第三判别函数在第三维中提供了最大的组分隔。以此类推。

1.2.2.1 原理
  1. 预测性LDA
    • 分组概率 同KNN,我们有

代替 ,则上式可以改写为:

LDA假定变量服从多元正态,且具有相等的协方差矩阵。我们用 代表合并的(pooled)组内样本协方差矩阵。 表示组 的样本均值。观测 的马氏平方距离是:

把上述样本估计值(2.1)带入多元正态密度函数,得到

将上式(2.2)代入 中,简化后得到:

通过上式(2.3)可以得到观测 属于组 的后验概率。

  • Cholesky分解与马氏变换

如果令 ,则

我们把(2.4)中的 称为马氏变换。

这样,我们可以通过两个步骤得到马氏距离:步骤一,利用马氏变换进行数据转换;步骤二,计算变换后数据的欧氏距离。

  1. 描述性LDA

假设有 个非零特征值,其中, 是组数, 是判别变量的个数。令 表示特征值,并且降序排列。令 表示相应的特征向量。

特征向量 乘以一个常数,仍旧是 特征值 的特征向量。特别地, 用作特征向量(判别函数),并称为非标准化规范判别函数,因为它们对应于没有标准化的判别变量。在vu向量的底部需要添加一个常数,这样在将vu向量应用于判别变量的线性组合时,可以保证所得变量的均值为零。

规范结构矩阵(canonical structure matrix)测量判别变量和判别函数之间的相关性,并在矩阵e(canstruct)中返回。规范结构矩阵等于 ,在第 行除以 ,其中 包含特征向量 作为列。

1.2.2.2 操作命令

语法结构:

discrim lda varlist [if] [in] [weight], group(groupvar) [options]

group(groupvar), 同knn,必选项,分组变量。

priors(priors), 设定各组先验概率

ties(ties),选择如何处理tie的情况。

notable,不报告分类表

lootable,显示留一法分类表。

1.2.2.3 估计后命令
estat anova,方差分析摘要
estat canontest,规范判别函数检验
estat classfunctions,分类函数
estat classtable,分类表
estat correlations,相关矩阵和P值
estat covariance,协方差矩阵
estat errorrate,错分率估计
estat grdistances,组均值的马氏平方距离和一般平方距离。
estat grmeans,组均值和转换后均值。
estat grsummarize,组统计摘要
estat list,分类表单
estat loadings,规范判别函数系数(载荷)
estat manova 多元方程分析表
estat structure,规范结构矩阵
estat summarize,估计样本摘要
loadingplot,画标准判别函数载荷
scoreplot, 画判别函数得分
screeplot,画特征值

estimates,估计结果目录
predict,预测分组结果和后验概率。

2. Stata操作实例

2.1 knn操作实例

  1. 生成训练数据 代码:
clear
set seed 98712321
set obs 500

generate x = 8*runiform() - 4
generate y = 8*runiform() - 4
generate group = 1
replace group = 2 if (y+2)^2 + (x+2)^2 <= 2
replace group = 2 if (y-2)^2 + (x-2)^2 <= 2
replace group = 3 if (y+2)^2 + (x-2)^2 <= 2
replace group = 3 if (y-2)^2 + (x+2)^2 <= 2

local rp : di %12.10f 2+sqrt(2)
local rm : di %12.10f 2-sqrt(2)

local functionplot ///
 (function y = sqrt(2-(x+2)^2) - 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = -sqrt(2-(x+2)^2) - 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = sqrt(2-(x-2)^2) + 2, lpat(solid) range(-`rm' `rp')) ///
 (function y = -sqrt(2-(x-2)^2) + 2, lpat(solid) range(-`rm' `rp')) ///
 (function y = sqrt(2-(x+2)^2) + 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = -sqrt(2-(x+2)^2) + 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = sqrt(2-(x-2)^2) - 2, lpat(solid) range( `rm' `rp')) ///
 (function y = -sqrt(2-(x-2)^2) - 2, lpat(solid) range( `rm' `rp'))

local graphopts ///
 aspectratio(1) legend(order(1 "group 1" 2 "group 2" 3 "group 3") rows(1))
 
 twoway (scatter y x if group==1) ///
 (scatter y x if group==2) ///
 (scatter y x if group==3) ///
 `functionplot' , `graphopts' name(original, replace) ///
 title("Training data")
图1 训练数据散点图
图1 训练数据散点图

上述代码生成3组数据,蓝色是1组,红色是2组,绿色是3组。由于1组围绕在2组和3组外围,不能使用线性判别方法,此处选择knn方法。

discrim knn x y, group(group) k(7) notable lootable priors(proportional)


predict cknn, looclass

local rp : di %12.10f 2+sqrt(2)
local rm : di %12.10f 2-sqrt(2)

local functionplot ///
 (function y = sqrt(2-(x+2)^2) - 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = -sqrt(2-(x+2)^2) - 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = sqrt(2-(x-2)^2) + 2, lpat(solid) range(-`rm' `rp')) ///
 (function y = -sqrt(2-(x-2)^2) + 2, lpat(solid) range(-`rm' `rp')) ///
 (function y = sqrt(2-(x+2)^2) + 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = -sqrt(2-(x+2)^2) + 2, lpat(solid) range(-`rp' -`rm')) ///
 (function y = sqrt(2-(x-2)^2) - 2, lpat(solid) range( `rm' `rp')) ///
 (function y = -sqrt(2-(x-2)^2) - 2, lpat(solid) range( `rm' `rp'))

local graphopts ///
 aspectratio(1) legend(order(1 "group 1" 2 "group 2" 3 "group 3") rows(1))
 
twoway (scatter y x if cknn==1 ) ///
 (scatter y x if cknn ==2) ///
 (scatter y x if cknn ==3) ///
 `functionplot', `graphopts' name(knn, replace) ///
 title("knn LOO classification")
图2 留一法生成的预测值
图2 留一法生成的预测值

从总错误率来看,knn的总错误率只有4%,取得了较为良好的分类效果。如果降低K值,总错误率可能会下降。当k=5时,knn的总错误率下降到3.4%;当k=3时,总错误率只有2%,但是留一法的错误率高达6.2%。也就是说,在k下降的同时,虽然普通的总错误率在下降,但是留一法的总错误率却上升了。这意味着,在k下降的时候,可能存在过拟合的现象。

2.2 lda操作实例

代码

clear mata
use https://www.stata-press.com/data/r16/twogroup, clear
discrim lda y x, group(group) notable
estat loadings, unstandardized
mat b = r(L_unstd)

mata
 c = st_matrix("b")
 b1 = -c[2,1]/c[1,1]
 b2 = -c[3,1]/c[1,1]
 b3 = c[1,1]/c[2,1]
 st_numscalar("b1",b1)
 st_numscalar("b2",b2)
 st_numscalar("b3",b3)
end

 twoway (sc y x if g == 1) (sc y x if g == 2) ///
 (function y = b1*x + b2,  range(5 65))(function y=b3*x-35,  range(40 80)), ///
 legend(order(1 "group 1" 2 "group 2" 3 "dividing line" 4 "projection line") ///
 size(*.7) rowgap(*.75) colgap(*.75) symxsize(*.75) symysize(*.75) colfirst) ///
 aspectratio(1) title(`"{fontface "楷体":Fisher线性判别分析}"') subtitle(`"{fontface "楷体":演示图}"')
图3 fisher判别
图3 fisher判别

从上图可以看到,fisher线性判别和knn有很大的不同。

参考文献

StataCorp LLC, Sturctural Equation Modeling Reference Manual[M]. A Stata Press Publication: Texas 2019

精益修身

2021/04/25  阅读:35  主题:嫩青

作者介绍

精益修身