机器学习－逻辑回归实例

上一讲中，我着重介绍了逻辑回归和最大似然估计的理论和推导，本文主要是在理论推导结果的基础上给出代码和实例。首先简要回顾一下上一讲中得到的两个重要结论：

似然函数为：\(\prod^m_{i=1} P(y^{(i)}=1｜x^{(i)}) \cdot P(y^{(i)}=0｜x^{(i)})\)；
逻辑回归梯度更新向量：\(\theta := \alpha \sum^m_{i=1} (g(x)-y) \cdot x\)。

本文主要从二项分布说起，直到构成逻辑回归的基本代价函数，以及目标函数的优化。

##1. 多重伯努利实验下的最大似然估计对于多重伯努利实验，我们可以获得产生产生序列的似然函数：\(\prod^n_{i=1} p^k \cdot (1-p)^{n-k}\)，对应的计算似然函数的代码如下：

likelihood <- function(sequence, p.parameter)
{
  likelihood <- 1
  num <- 0
  for (i in 1:length(sequence))
  {
    if (sequence[i] == 1)
    {
      #注意这里没有求对数，会产生问题，接下来也会解释这些问题
      likelihood <- likelihood * p.parameter
      num <- num + 1
    }
    else
    {
      likelihood <- likelihood * (1 - p.parameter)
    }
  }
  likelihood <- likelihood
  return(likelihood)
}

假使多重伯努利实验的概率\(p\)为0.8，我们可以采用R中的rbinom函数生成对应多重伯努利实验序列，很明显，产生该序列的最优\(p\)值为序列中出现1个数的比例。

#生成原始伯努利序列，最优值明显是序列的均值
p.parameter <- 0.8
sequence <- rbinom(10, 1, p.parameter)
#计算从0到1的概率值对应的序列的似然函数值
possible.p <- seq(0, 1, by = 0.01)
VectorP <- sapply(possible.p, function (p) {likelihood(sequence, p)})
##绘制似然函数图，打印最大的似然函数值，即产生序列的最可能概率
plot(possible.p, VectorP, type = "l", lty = 1, lwd = 1)
print(sprintf("产生的序列为%g", sequence))
print(sprintf("最优值为序列均值：%g", possible.p[which(VectorP == max(VectorP))]))

通过绘制的似然函数曲线或者我们的直观印象可以猜到产生该序列的最优p值为0.9。接下来，我们采用最大似然估计对产生序列的最大可能的似然函数的概率进行求解，在R中对于一维变量，可以选用optimize函数来逼近最优值，运行下述代码，可以获得优化结果为0.899998，结果和我们常识认为的结果一致：

#采用optimize优化方法寻优
mle.results <- optimize(function(p) {likelihood(sequence, p)},
                        interval = c(0, 1),
                        maximum = TRUE)
print(sprintf("optimize()函数优化结果为：%g", mle.results$maximum))

另外，我们通过机器学习的学习可以了解到，优化算法以及最大似然估计会随着数据量的增加而使得估计结果更为准确，因此，接下来我们来计算在不同数据量下估计误差，具体实现代码及结果如下：

error.behavior <- data.frame()
findSample <- seq(100, 2500, 100)
#计算数据量从100到2500大小的MLE估计误差
for (n in findSample)
{
  sequence <- rbinom(n, 1, p.parameter)
  likelihood.results <- optimize(function(p) {likelihood(sequence, p)},
                                 interval = c(0, 1),
                                 maximum = TRUE)
  true.mle <- mean(sequence)
  likelihood.error <- true.mle - likelihood.results$maximum
  error.behavior <- rbind(error.behavior,
                          data.frame(N = n,
                                     Error = likelihood.error,
                                     Algorithm = 'Likelihood'))
}
#绘制误差曲线
plot(error.behavior[ ,1], error.behavior[ ,2], type = "l",
     xlab = "number of sample", ylab = "estimation error", 
     main = "Error Behavior of Naive MLE")
print(sprintf("the sample is %g, and the likelihood is %g", 
              2500, likelihood.results$maximum))

从上述朴素MLE实验结果发现当样本量大于1500时，估计误差陡然增加，该结果与优化思想中的渐进逼近理论相悖，应该是代码中存在错误。经发现，当样本量为2500时，似然值为0.9999，接近于1，分析其原因发现，应当是浮点型小数连乘后造成数特别小接近于0，导致计算机计算精度不够，结果为NaN。理想的改进方案是求对数！！这是不是也从代码层面上解释了为什么我们要把优化函数转换为凸函数，为什么似然函数要求对数似然函数呢？？所以很明显，我们需要对似然函数的计算方式进行改变，具体代码如下：

#要求对数似然函数
log.likelihood <- function(sequence, p)
{
  log.likelihood <- 0
  for (i in 1:length(sequence))
  {
    if (sequence[i] == 1)
    {
      log.likelihood <- log.likelihood + log(p)
    }
    else
    {
      log.likelihood <- log.likelihood + log(1 - p)
    }
  }
  return(log.likelihood)
}

#重新估计不同样本量下的误差
sequence <- rbinom(10, 1, p.parameter)
#计算从0到1的概率值对应的序列的似然函数值
possible.p <- seq(0, 1, by = 0.001)
VectorLogP <- sapply(possible.p, function (p) {log.likelihood(sequence, p)})
##绘制似然函数图，打印最大的似然函数值，即产生序列的最可能概率
plot(possible.p, VectorLogP, type = "l", lty = 1, lwd = 1,
     xlab = "probability", ylab = "log-likelihood", main = "Log Likelihood of function")

log.error.behavior <- data.frame()
findLogSample <- seq(100, 2500, 100)
for (n in findLogSample)
{
  sequence <- rbinom(n, 1, p.parameter)
  log.likelihood.results <- optimize(function(p) {log.likelihood(sequence, p)},
                                 interval = c(0, 1),
                                 maximum = TRUE)
  true.mle <- mean(sequence)
  log.likelihood.error <- true.mle - log.likelihood.results$maximum
  log.error.behavior <- rbind(log.error.behavior,
                          data.frame(N = n,
                                     Error = log.likelihood.error,
                                     Algorithm = 'Likelihood'))
}
plot(log.error.behavior[ ,1], log.error.behavior[ ,2], type = "l",
     xlab = "number of sample", ylab = "estimation error", 
     main = "Error Behavior of Logriathm MLE")
print(sprintf("the sample is %g, and the log.likelihood is %g", 
              2500, likelihood.results$maximum))

从上面结果可以看出，预测误差基本维持在0附近，所以通过求对数会使得MLE估计更为准确。

##2. 最大似然估计算法伪代码通过上面多重伯努利实验下的最大似然估计的实例分析，可以总结出最大似然估计的基本伪代码，如下：

#最大似然函数的伪代码
log.likelihood <- function(sequence.as.data.frame, likelihood.function, parameters)
{
  log.likelihood <- 0
  
  for (i in 1:nrow(sequence.as.data.frame))
  {
    log.likelihood <- log.likelihood + log(likelihood.function(
      sequence.as.data.frame[i,], parameters))
  }
  
  return(log.likelihood)
}

##3. 最大似然估计求解logistic回归对于logistic回归而言，分类结果的表达式如下： \(p(y=1|x,\theta)=exp(\theta∗x)/(1+exp(\theta∗x))\)和\(p(y=0|x,\theta)=1/(1+exp(\theta∗x))\)。所以，为采用MLE求解logistics回归参数，需要定义似然函数。为获得产生样本的最大可能，根据贝叶斯误差，对于所有样本，\(p(y=1|x,\theta)\)或者\(p(y=0|x,\theta)\)要尽可能最大。所以对于单一样本，似然函数可定义为:

\[ll=−y∗p(y=1|x,\theta)−(1−y)∗p(y=0|x,\theta)\]

之前博客中提到的数值分析中的optimize()函数适用于一维的优化问题，多维度优化问题应当采用optim()函数。本实例采用iris数据做二分类，并分别在给定梯度和未给定梯度的情况下，选用两种优化方法的求解过程。具体代码实例如下：

#更改iris数据类标签，改为二分类问题
irisData <- iris
irisData$Species <- factor(irisData$Species, labels = c('0', '1', '2'))
irisData$Species <- as.numeric(irisData$Species)
irisData[which(irisData$Species != 1), ]$Species = 0

#logistic回归的似然函数
lr.log.likelihood <- function(x, para){
  log.likelihood <- 0
  RlenX <- length(x[1,])
  intervalX <- rep(1, RlenX)
  x <- cbind(intervalX, x)
  ClenX <- length(x[1,])
  
  for (i in 1:length(x[ ,1]))
  {
    log.likelihood <- log.likelihood - 
      x[i,ClenX] * log(exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para)) / (1+exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para)))) - (1-x[i,ClenX]) * log(1/(1+exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para))))
  }  
  log.likelihood <- log.likelihood / length(x[,1])
  return(log.likelihood)
}

#选用单纯形算法求解最优值
##给定初始参数
initPara <- c(0.5, 3, 1.0, 2.0, 1)
lr.mle.results <- optim(par = initPara,fn = function(p) {lr.log.likelihood(irisData, p)}, 
                        method = "Nelder-Mead",  #或者采用SANN方法
                        control = list(trace = TRUE, maxit = 100))

#预测类别标签
final.mle.para <- lr.mle.results$par
mle.predict <- matrix(0 ,length(irisData$Species))
Rlen <- length(irisData[1,])
interval.data <- rep(1, Rlen)
bindMleData <- cbind(interval.data, irisData)
beta.mle.data <- as.matrix(bindMleData[,1:(length(bindMleData[1,])-1)]) %*% as.matrix(final.mle.para)
mle.probs <- exp(beta.mle.data) / (1 + exp(beta.mle.data))
mle.predict[mle.probs > 0.5] = 1
mle.predict[mle.probs <= 0.5] = 0
print(table(observed=irisData$Species, predicted=mle.predict))

通过优化目标函数（对数似然函数），并构建分类logistic回归模型，预测结果准确率为100%。接下来，因为我们可以求得逻辑回归代价函数的梯度，所以也可以选用bfgs算法求解极值。具体实现代码如下，通过梯度下降算法，我们仍可以获得100%的分类准确率：

#给定梯度，采用BFGS算法求解
bfgs.likelihood <- function(x, para){
  log.likelihood <- 0
  RlenX <- length(x[1,])
  intervalX <- rep(1, RlenX)
  x <- cbind(intervalX, x)
  ClenX <- length(x[1,])
  
  for (i in 1:length(x[ ,1]))
  {
    log.likelihood <- log.likelihood - 
      x[i,ClenX] * log(exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para)) / (1+exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para)))) - (1-x[i,ClenX]) * log(1/(1+exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para))))
  }  
  log.likelihood <- log.likelihood / length(x[,1])
  return(log.likelihood)
}

gr.log.likelihood <- function(x, para){
  RlenX <- length(x[1,])
#  p <- 1/(1 + exp(as.matrix(-x) %*% as.matrix(para)))
#  -(crossprod(as.matrix(x),as.matrix(x[,RlenX] - p)))
   gr.likelihood <- rep(0, length(para))
   intervalX <- rep(1, RlenX)
   x <- cbind(intervalX, x)
   ClenX <- length(x[1,])
   for (j in 1:length(para)){
     for (i in 1:length(x[ ,1]))
     {
       gr.likelihood[j] <- gr.likelihood[j] + 
         x[i,j] * (exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para))/(1 + exp(as.matrix(x[i,1:(ClenX-1)]) %*% as.matrix(para))) - x[i,ClenX])
     }
   }
   gr.likelihood <- gr.likelihood / length(x[,1])
#   return(gr.likelihood)
}

lr.bfgs.results <- optim(par = initPara,
                        fn = function(p) {bfgs.likelihood(irisData, p)},
                        gr = function(p) {gr.log.likelihood(irisData, p)},
                        method = "L-BFGS-B", #L-BFGS-B
                        control = list(trace = TRUE, maxit = 100))

##4. 总结通过上述实例的分析，我们可以从代码层面剖析最大似然估计以及逻辑回归的问题，并证明了最大似然估计需要选用对数似然变换。

机器学习－逻辑回归实例

预计阅读时间： 57分钟

总计阅读次数：次