本次为学生时期所写的实验报告，代码程序为课堂学习和自学，对网络程序有所参考，如有雷同，望指出出处，谢谢！

基础知识来自教材：李航的《统计学习方法》

本人小白，仍在不断学习中，有错误的地方恳请大佬指出，谢谢！

本次实验将分别采用梯度下降法和牛顿法进行二项逻辑斯蒂回归模型参数估计

一、代码实现（R语言）

1.梯度下降法

#预备代码：
#计算似然函数：对i求和( yi*(w^T*x)-log(1+exp(w^T*x)) )
likelihood_fuc<-function(y,X,w){N=nrow(X)save=rep(NA,N)for (i in 1:N) {save[i]=y[i]*t(w)%*%X[i,]-log(1+exp(t(w)%*%X[i,]))}sum=sum(save)sum
}#梯度下降法
my_gradiant_descent<-function(X,y,w,ddmax,limit,lambda){  #X为n*p+1矩阵(n次观测，p个特征，最后一列全为1)#y为n维列向量#w输入初值w0(w为p+1维向量)#ddmax为最大迭代次数#limit为给定阈值#lambda为一系列步长的预设值    N=dim(X)[1]    #N为样本数pai=rep(NA,N) #π(x)为n维列向量，记录每个给定xi（xi为p维向量）下y=1的概率i=1while (i<=ddmax){#给定x、w时y=1的概率π（x）
pai=t(exp(t(w)%*%t(X))/(1+exp(t(w)%*%t(X)))) #pai为n×1的列向量#计算此时梯度X^T*(π-y)grad=t(X)%*%(pai-y)grad_norm=sqrt(sum(grad^2)) #梯度的模长#迭代的判断条件及其处理if(grad_norm>=limit){#选取最适宜步长（寻找给定w下使似然函数最大的步长）result=rep(NA,length(lambda))for (k in 1:length(lambda)) {result[k]=likelihood_fuc(y,X,w-lambda[k]*grad)}lambdabest=lambda[which.max(result)]#更新ww=w-lambdabest*grad#进入下一次迭代i=i+1}else{break    #若不满足阈值要求，则退出循环}}w
}

2.牛顿法

my_Newton_method<-function(X,y,w,ddmax,limit){  #X为n*p+1矩阵(n次观测，p个特征，最后一列全为1)#y为n维列向量#w输入初值w0(w为p+1维向量)#ddmax为最大迭代次数#limit为给定阈值N=dim(X)[1]  #N为样本数pai=rep(NA,N)
#π(x)为n维列向量，记录每个给定xi（xi为p维向量）下y=1的概率D=matrix(0,ncol=N,nrow = N)    #D为n×n的空矩阵，用于后续海塞矩阵的便捷计算 i=1while (i<=ddmax){#给定x、w时y=1的概率π（x）for(j in 1:N){pai[j]=exp(t(w)%*%X[j,])/(1+exp(t(w)%*%X[j,]))D[j,j]=pai[j]*(1-pai[j])
#D为对角矩阵，第(i,i)元素为π(xi)*(1-π(xi));其中xi为每一次观测的值，为p+1维向量}#计算此时梯度grad=t(X)%*%(pai-y)grad_norm=sqrt(sum(grad^2)) #梯度的模长#迭代的判断条件及其处理if(grad_norm>=limit){#海塞矩阵求逆invH=solve(t(X)%*%D%*%X)#更新ww=w-invH%*%grad#进入下一次迭代i=i+1}else{break}}w
}

3. S-折交叉验证函数编写

###梯度下降法交叉验证
Stest_gradiant_descent<- function(X,y,S,ddmax,limit,lambda,up,down){    #X为n*p+1矩阵(n次观测，p个特征，最后一列全为1)#y为n维列向量#S为折数#ddmax为最大迭代次数；limit为阈值#up、down为预设阈值：即概率P(Y=1|x)大于up时认定y为1，P(Y=1|x)小于给定down时认定y为0#数据的行数n = dim(X)[1]   percentage = rep(NA,S)#进行1到S的循环，总长度为n（如1,2,3,4,1,2,3,4,1,2...），并进行打乱CV.ID = sample(rep(1:S,length.out = n))  for (j in 1:S){
X.test = X[CV.ID==j,]
#相当于从被分为S份的数据中随机抽取1份作为测试数据集，循环S次y.test = y[CV.ID==j]
X.train = X[CV.ID!=j,]
#相当于从被分为S份的数据中随机抽取（S-1）份作为训练数据集，循环S次y.train = y[CV.ID!=j]w.hat=my_gradiant_descent(X.train,y.train,w=c(0,0,0,0,0),
ddmax,limit,lambda)                                    #测试数据集中y的预测值y.hat.prop=exp(t(w.hat)%*%t(X.test))/
(1+exp(t(w.hat)%*%t(X.test)))
#y.hat.prop为n维行向量，每个元素代表Y=1的概率y.hat.prop[y.hat.prop>=up]=1    #预设概率大于up时为1，概率小于down时为0y.hat.prop[y.hat.prop<down]=0y.hat=y.hat.prop#由准确率来判断误差大小percentage[j] = sum(y.hat==y.test)/length(y.test)}correct = mean(percentage)error = 1-mean(percentage)plot(1:S,percentage,type="l")   #绘制准确率的折线图list(correct = correct,error = error,percentage = percentage)}###牛顿法交叉验证
Stest_Newton<- function(X,y,S,ddmax,limit,up,down){    #X为n*p+1矩阵(n次观测，p个特征，最后一列全为1)#y为n维列向量#S为折数#ddmax为最大迭代次数；limit为阈值#up、down为预设阈值：即概率P(Y=1|x)大于up时认定y为1，P(Y=1|x)小于给定down时认定y为0#数据的行数n = dim(X)[1]   percentage = rep(NA,S)#进行1到S的循环，总长度为n（如1,2,3,4,1,2,3,4,1,2...），并进行打乱CV.ID = sample(rep(1:S,length.out = n))  for (j in 1:S){
X.test = X[CV.ID==j,]
#相当于从被分为S份的数据中随机抽取1份作为测试数据集，循环S次y.test = y[CV.ID==j]
X.train = X[CV.ID!=j,]
#相当于从被分为S份的数据中随机抽取（S-1）份作为训练数据集，循环S次y.train = y[CV.ID!=j]w.hat=my_Newton_method(X.train,y.train,w=c(0,0,0,0,0),
ddmax,limit)   #测试数据集中y的预测值y.hat.prop=exp(t(w.hat)%*%t(X.test))/
(1+exp(t(w.hat)%*%t(X.test)))
#y.hat.prop为n维行向量，每个元素代表Y=1的概率y.hat.prop[y.hat.prop>=up]=1    #预设概率大于up时为1，概率小于down时为0y.hat.prop[y.hat.prop<down]=0y.hat=y.hat.prop#由准确率来判断误差大小percentage[j] = sum(y.hat==y.test)/length(y.test)}correct = mean(percentage)error = 1-mean(percentage)plot(1:S,percentage,type="l") #绘制准确率的折线图list(correct = correct,error = error,percentage = percentage)}

二、检验

1.输入

##########数据处理####
##选择数据：iris数据集
iris <- read.csv("F:/大三下/数据挖掘/UCI数据集/Iris/iris.data",header = FALSE)
y<-iris[,5]#由于为二项逻辑回归，考虑只选择y的种类为"Iris-setosa"和"Iris-virginica"的数据集
split=split(iris,factor(y),drop = TRUE)
#使用split函数将数据按特征列y的各取值进行分割成多个列表
#矩阵X
X=rbind(split[[1]],split[[3]])[,-5]
rownames(X)=NULL
X=cbind(X,1)
X=as.matrix(X)
X#种类列y
y=rbind(split[[1]],split[[3]])[5]
rownames(y)=NULL
y[y=="Iris-setosa"]=1
y[y=="Iris-virginica"]=0
#将y两个取值"Iris-setosa"、"Iris-virginica"分别用1、0表示
y=as.numeric(as.matrix(y))
y#此时X与y的类型均为矩阵，X的规格为n×(p+1)
################检验过程####
#给定概率精度0.9、0.1
Stest_gradiant_descent(X,y,S=5,ddmax=100,limit = 0.01,
lambda=seq(0.001,0.005,0.001),
up=0.9,down=0.1)Stest_Newton(X,y,S=5,ddmax=100,limit = 0.01,up=0.9,down=0.1)#给定概率精度0.9999、0.1
Stest_gradiant_descent(X,y,S=5,ddmax=100,limit = 0.01,
lambda=seq(0.001,0.005,0.001),
up=0.9999,down=0.1)Stest_Newton(X,y,S=5,ddmax=100,limit = 0.01,up=0.9999,down=0.1)

2.输出结果与分析

（1）给定概率精度0.9、0.1

①梯度下降法

>Stest_gradiant_descent(X,y,S=5,ddmax=100,limit = 0.01,lambda=seq(0.001,0.005,0.001),up=0.9,down=0.1)

$correct

[1] 1

$error

[1] 0

$percentage

[1] 1 1 1 1 1

分析：由结果可知：给定概率精度0.9、0.1下，使用梯度下降法估计的w进行预测时，平均预测准确率为100%，五次交叉验证的预测准确率都为100%，预测效果很好

②牛顿法

> Stest_Newton(X,y,S=5,ddmax=100,limit = 0.01,up=0.9,down=0.1)

$correct

[1] 1

$error

[1] 0

$percentage

[1] 1 1 1 1 1

分析：由结果可知：给定概率精度0.9、0.1下，使用牛顿法估计的w进行预测时，平均预测准确率为100%，五次交叉验证的预测准确率都为100%，预测效果很好。

（2）给定概率精度0.9999、0.1

①梯度下降法

$correct

[1] 0.5

$error

[1] 0.5

$percentage

[1] 0.65 0.40 0.50 0.50 0.45

分析：由结果可知：给定概率精度0.9999、0.1下，使用梯度下降法估计的w进行预测时，平均预测准确率为50%，五次交叉验证的预测准确率均徘徊在50%之间，预测效果变差。

②牛顿法

> Stest_Newton(X,y,S=5,ddmax=100,limit = 0.01,up=0.9999,down=0.1)

$correct

[1] 0.95

$error

[1] 0.05

$percentage

[1] 0.95 0.90 0.95 1.00 0.95

分析：由结果可知：给定概率精度0.9999、0.1下，使用牛顿法估计的w进行预测时，平均预测准确率为95%，五次交叉验证的预测准确率均在90%及以上，预测效果仍然较好。

三、小结

本次实验主要编写了梯度下降法和牛顿法对逻辑斯地回归似然函数参数进行估计。首先需熟练掌握二项逻辑斯地回归的似然函数的建立；其次在应用算法时需明确似然函数对参数w的一阶梯度和二阶梯度。梯度下降法牛顿法都是迭代算法，通过选取适当的初值x(0)，不断迭代，更新x的值，进行目标函数的最小化，直至收敛。负梯度方向是使函数下降最快的方向，所以在迭代的每一步，以负梯度方向更新x的值从而达到减少函数值的目的。而梯度下降法与牛顿法的不同之处在于：梯度下降法使用一阶泰勒逼近；而牛顿法使用了二阶泰勒逼近。这使得两种算法对x的更新方法不同。在正式的函数编写和交叉验证中，可以发现牛顿法的运行速度较梯度下降法快，在适当精度下，通过牛顿法与梯度下降法估计的参数进行预测的准确率都很好，但在高精度情况下，梯度下降法估计的参数预测准确率有明显下降，而牛顿法则维持原来的较高的准确率水平。

数据挖掘-二项逻辑斯蒂回归模型算法的R实现相关推荐

逻辑斯蒂分布模型、二项逻辑斯蒂回归模型、多项逻辑斯蒂回归模型
一.逻辑斯蒂分布/回归模型模型描述的是一种什么样的事件或现象: 设X是连续随机变量,X服从逻辑斯蒂回归分布是指X具有下列分布函数和密度函数: 附上逻辑斯蒂分布的密度函数与分布函数,如下: 物理含义, ...
回归分析（三）二项逻辑斯蒂回归模型
回归分析(三)二项逻辑斯蒂回归学了一段时间突然又遇到逻辑斯蒂回归,结果发现已经忘完了,所以今天重新梳理一下. (1)逻辑斯蒂分布先看一下逻辑斯蒂分布函数F(x)F(x)F(x),其概率密度 ...
逻辑斯蒂回归模型——逻辑斯蒂分布、二项逻辑斯蒂回归模型、参数估计与多项逻辑斯蒂回归
本笔记整理自李航老师<统计学习方法>第二版第六章逻辑斯蒂回归是统计学习中经典的分类方法. 逻辑斯蒂分布 F(x)=P(X≤x)=11+e−(x−μ)/γF(x) = P(X\leq x ...
二项逻辑斯蒂回归（逻辑回归）
二项逻辑斯蒂回归 1 定义二项逻辑斯蒂回归模型是如下的条件概率分布: P(Y=1∣x)=exp⁡(ω⋅x+b)1+exp⁡(ω⋅x+b)P(Y=0∣x)=11+exp⁡(ω⋅x+b)P(Y=1 | ...
用二项逻辑斯蒂回归解决二分类问题
逻辑斯蒂回归: 逻辑斯蒂回归是统计学习中的经典分类方法,属于对数线性模型.logistic回归的因变量可以是二分类的, 也可以是多分类的基本原理 logistic 分布折X是连续的随机变量,X服从 ...
理解逻辑斯蒂回归模型
逻辑斯蒂回归是一个非常经典的二项分类模型,也可以扩展为多项分类模型.其在应用于分类时的过程一般如下,对于给定的数据集,首先根据训练样本点学习到参数w,b;再对预测点分别计算两类的条件概率,将预测点判为 ...
回归分析：逻辑斯蒂回归模型，可视化分类决策边界
文章目录逻辑斯蒂回归模型逻辑斯蒂回归模型python案例逻辑斯蒂回归模型前面的例子都是在用线性模型解决回归任务,那么线性模型能否完成分类任务呢?相较于回归任务,分类任务的预测值是离散的,比如二 ...
回归分析：逻辑斯蒂回归模型，非线性分类任务案例
逻辑斯蒂回归模型上节内容,我们得到逻辑斯蒂回归模型,需要注意虽然名称叫做回归模型,但是它是解决分类任务,公式如下: y=11+e−(wTX+b)y = \frac {1}{1+e^{-(w^T X ...
逻辑斯蒂回归分类算法
逻辑斯蒂回归分类算法首先来看一个线性回归来进行分类的问题: 怎样判断肿瘤是否恶性? 很明显线性回归用于分类问题无法处理边界点的位置. 同时,线性回归健壮性不够,一旦有噪声,立刻"投降&qu ...

数据挖掘-二项逻辑斯蒂回归模型算法的R实现

一、代码实现（R语言）

1.梯度下降法

2.牛顿法

3. S-折交叉验证函数编写

二、检验

1.输入

2.输出结果与分析

三、小结

数据挖掘-二项逻辑斯蒂回归模型算法的R实现相关推荐

最新文章

热门文章