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_data
plot_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_data
traindata$x0 <- 1
traindata <- 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个
#Visualize
classify_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.765
Hession_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个
#Visualize
classify_visualization(plot_data,0.62,7.38,"牛顿法求解")

结果错分四个