Logistic回归

Logistic回归

Part1:优缺点

优：直接对分类的可能性进行建模，无需事先假设数据分布。对率函数是任意阶可导的凸函数，方便求解 缺：对自变量的多重共线性比较敏感，预测结果呈S型分布，变化率越来越小，导致很多区间的变量变化对目标概率没有区分度 优点：计算代价不高，易于理解和实现。 缺点：容易欠拟合，分类精度可能不高。 适用数据类型：数值型和标称型数据。 from 周志华《机器学习》 1、无需事先假设数据分布，避免假设分布不准确所带来的问题 2、可得到近似概率预测，可以辅助决策 3、对率函数是任意阶可导的凸函数，有很好的数学性质，可以直接用于求解最优值[貌似是损失函数是凸函数，是书里写错了吗]

Part2：理论推导

1、Logistic Function

其中：

2、回归模型

注意最后一个合并的形式，方便求解~

3、最大似然估计

还可以进一步化简

4、目标函数为凸函数

（1）凸函数的定义：

（2）凸函数等价条件

（3）凸函数的运算

（4）目标函数

前面的极大似然估计需要求最大值。还是从常用的损失函数考虑，逻辑上更清晰一点。

考虑两种错判的损失：

如果样本实际是1，则返回的h(x)

损失函数：

损失函数的凹凸性：

由凸函数的性质可知，两个凸函数的和是凸函数。那么对于损失函数而言，我们只需要分析单个样本损失函数的凹凸性即可.

可见都是凸函数，那么目标函数存在唯一极小值点。

画个图感受下损失函数：

opar <- par(no.readonly = TRUE)par(mfrow=c(1,2))curve(expr = -log(1/(1+exp(-x))),from = -10,to=10,n=1000,xlab="z of y=1",ylab="cost")curve(expr = -log(1-1/(1+exp(-x))),from = -10,to=10,n=1000,xlab="z of y=0",ylab="cost")par(opar)

5、参数求解

（1）R自带glm函数求解：

# 机器学习实战数据 ----------------------------------------------------------------logistic_data <- read.table("clipboard")colnames(logistic_data) <- c("x1","x2","y")##预测为versicolor的概率glm.fit <- glm(y~x1+x2, data=logistic_data,family="binomial")summary(glm.fit)pro <- predict(glm.fit,logistic_data[,1:2],type = "response")pre_label <- ifelse(pro>0.5,1,0)table(pre_label,logistic_data$y)#错判了5个样本plot_data <- logistic_dataplot_data$y <- as.factor(plot_data$y)#传入参数:数据集、斜率、截距、标题classify_visualization <- function(plot_data,slope,intercept,title){ library(ggplot2) p <- ggplot(data=plot_data,aes(x=x1,y=x2,color=y,shape=y))+ geom_point(size=2.5)+ scale_colour_manual(values = c( "red", "blue"))+#指定点的颜色 scale_shape_manual(values = c( 16, 15))+#指定点的形状16圆、15方形 geom_abline(slope =slope ,intercept = intercept)+#画分割超平面 labs(title = title)+#add title theme( plot.title = element_text( size = 15, hjust =0.5, vjust = 1, color="black" ) #改变标题的位置、颜色、字体大小 ) return(p)}classify_visualization(plot_data,0.62,7.38,"R自带glm函数求解")

结果：错判5个样本

（2）随机梯度上升法

求导自己个算

# 随机梯度梯度上升求参 ------------------------------------------------------------------###没有用矩阵乘法，代码写的冗长。。。#传入特征、参数logistic_f <- function(x, theta) { z <- exp(sum(x * theta)) return(z / (1 + z))}# 只针对这份数据做计算上面的认证# 传入参数:trainset训练的数据,数据框；theta参数Likelihood <- function(traindata,Theta){ a <- traindata[1]*Theta[1]+traindata[2]*Theta[2]+traindata[3]*Theta[3] b <- exp(a)/(1+exp(a)) c <- b^traindata[4]*(1-b)^(1-traindata[4]) L <- cumprod(c)[length(c)] return(L)}## 随机梯度上升求解traindata <-logistic_datatraindata$x0 <- 1traindata <- traindata[,c(4,1:3)]J_r <- c()#似然值Theta <- c(1, 1, 1)#参数初始值Alpha <- 0.01#步长for(i in 1:500) { sample_seed <- sample(1:nrow(traindata), 1, replace = F) trainset <- as.matrix(traindata[sample_seed,])[1,] y <- trainset[4] Theta <- Theta+Alpha*(y-logistic_f(trainset[1:3],Theta))*trainset[1:3] print(Theta) L <- Likelihood(traindata,Theta) print(L) J_r <- c(J_r,L)}#后验概率x <- as.matrix(traindata)[,-4]theta <- as.matrix(Theta)z <- exp(x %*% theta)zz <- (z / (1 + z))table(traindata$y,ifelse(zz>0.5,1,0))#错分7个样本，但是随机梯度的结果变动太大了#可视化lassify_visualization(plot_data,1.32,5.86,"随机梯度上升求解")

结果：错判7个样本

（3）批梯度上升法

涉及到矩阵求导，还是要细心点

# 批梯度上升 -------------------------------------------------------------------#log_f函数传入两个参数：观测矩阵、系数矩阵log_f <- function(x, theta) { z <- exp(x %*% theta) return(z / (1 + z))}Stochastic_gradient_descent <- function(x,y){ x <- as.matrix(x) y <- as.matrix(y) x_transpose <- t(x) Alpha <- 0.001 Theta <- matrix(c(1, 1, 1),byrow = F) for(i in 1:500){ error <- y-log_f(x,Theta) Theta <- Theta+Alpha*x_transpose%*%error } return(Theta)}x=traindata[,1:3]y=traindata[,4]theta <- Stochastic_gradient_descent(x,y)#后验概率x <- as.matrix(traindata)[,-4]z <- exp(x %*% theta)zz <- (z / (1 + z))table(traindata$y,ifelse(zz>0.5,1,0))#错分4个#Visualizeclassify_visualization(plot_data,0.77,6.65,"批梯度上升求解")

结果：错分4个

（4）牛顿法

理论还是不太理解，感觉没有梯度下降明显啊~这脑子哟

两个问题：

1、二阶泰勒展开用”=”，并且对等号两边求导，合理不？

2、是否是假设Hesse matrix正定？否则一阶导为0，不是充要条件啊。凸优化的教材上说只有在目标函数为凸函数时候，该条件才是充要的？？？

自问自答：logloss是高阶可导连续凸函数~至于证明，台大的视频课里似乎有，看后再补吧~

hession矩阵的求解要了老命。。。

R代码，结果和R一毛一样~

# 牛顿法 ---------------------------------------------------------------------#log_f函数传入两个参数：观测矩阵、系数矩阵log_f <- function(x, theta) { z <- exp(x %*% theta) return(z / (1 + z))}x=traindata[,1:3]#特征矩阵y=traindata[,4]#因变量x_square <- as.matrix(x)*as.matrix(x)# theta <- matrix(1:3)# x_feature <- as.matrix(x)#R软件绝对误差和5.765Hession_Matrix <- function(theta,x_square,x_feature){ y1 <- -x_feature%*%theta y2 <- exp(y1) y <- -y2/(1+y2)^2 h11 <- t(x_square[,1])%*%y h22 <- t(x_square[,2])%*%y h33 <- t(x_square[,3])%*%y h12 <- t(x_feature[,1]*x_feature[,2])%*%y h13 <- t(x_feature[,1]*x_feature[,3])%*%y h23 <- t(x_feature[,2]*x_feature[,3])%*%y hession_matrix <- matrix(c(h11,h12,h13,h12,h22,h23,h13,h23,h33),3,3) return(hession_matrix)}Newton_descent <- function(x,y){ x_feature <- as.matrix(x) y <- as.matrix(y) x_feature_transpose <- t(x_feature) Alpha <- 0.5 Theta <- matrix(c(0, 0, 0),byrow = F) for(i in 1:25){ error <- y-log_f(x_feature,Theta) print(sum(abs(error))) gradient <- x_feature_transpose%*%error print(gradient) hession_inverse <- solve(Hession_Matrix(Theta,x_square,x_feature)) Theta <- Theta-Alpha*hession_inverse%*%gradient print(Theta) } return(Theta)}theta <- Newton_descent(x=traindata[,1:3],y=traindata[,4])#后验概率x <- as.matrix(traindata)[,-4]z <- exp(x %*% theta)zz <- (z / (1 + z))table(traindata$y,ifelse(zz>0.5,1,0))#错分5个#Visualizeclassify_visualization(plot_data,0.62,7.38,"牛顿法求解")

结果错分四个

版权声明：本文内容由网络用户投稿，版权归原作者所有，本站不拥有其著作权，亦不承担相应法律责任。如果您发现本站中有涉嫌抄袭或描述失实的内容，请联系我们jiasou666@gmail.com 处理，核实后本网站将在24小时内删除侵权内容。