大白话介绍逻辑回归牛顿法求解+python实现

首先，对于“受多个因素影响会产生两种结果“的这类问题，我们想要训练一种模型，通过这种模型能够实现任意输入一个n维向量，可以输出该类问题是哪一类。

对于上述的情况，我们可以考虑采用逻辑回归的思想来解决。

逻辑回归是输入一个n维向量，通过加权求和，经过激活函数的“打包”，综合得出一个“概率值”，根据该“概率值”是否达到某个阈值来判断结果是1还是0.

具体来说：（假如我的自变量是5维的）

 注意：如果你的自变量是n维的，那么你要输入的其实是n+1维的，其中有个偏置不可以忽略

我们注意到，要想解决这个问题，我们目前未知的是这个n+1维的theta向量，因此如何求得这个theta向量是至关重要的。

求theta向量这里要引入一个极大似然的概念，具体来说：

任意一个“概率”的公式为：

则它的似然函数为：

经简化：

这里我们要求解似然函数的最大值，进一步问题转化为最优化问题，又进一步转化为其导函数求根的问题——牛顿法！

牛顿法是一种通过迭代近似求根的方法，使用函数f (x)的泰勒级数的前面几项来寻找方程f (x) = 0的根。用牛顿迭代法解非线性方程，是把非线性方程f(x) = 0线性化的一种近似方法。

把f(x)在点x0的某邻域内展开成泰勒级数：

取其前两项可以得到：

进而可以得到牛顿法迭代关系式：

注意：这里的x是我们的theta向量，这里的分母虽然在泰勒里是一阶导，但事实上是两阶导，因为f(x)本身就是一阶导。

下面是python源码，亲写亲测可运行无报错：

#利用牛顿法求解逻辑回归，然后画图（详细版）
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
data=pd.read_csv("E:\\caffe\\study\\2_data.csv",sep=",")
#整理好数据
train_x=data.iloc[0:100,0:2]
train_y=data.iloc[0:100,2:3]
theta=np.mat([1,1,1])#初始theta
times=30#定义迭代次数
#定义函数
def h(theta,x):#LR假设函数
    return(1.0/(1+np.exp(-1*(np.dot(theta,x)))))
def j(train_x,train_y,theta):#一阶导函数
    a=np.mat([0,0,0])
    # train_x=pd.concat([DataFrame([1 for i in range(len(train_x))]),train_x],1)
    train_x["x3"]=[1.0 for i in range(len(train_x))]
    for i in range(len(train_x)):
        # a = a - (h(theta, np.mat(train_x.iloc[i, :]).T) - train_y.iloc[i, 0]) * np.mat(train_x.iloc[i, :])
        a=a+ (h(theta,np.mat(train_x.iloc[i,:]).T)-train_y.iloc[i,0])*np.mat(train_x.iloc[i,:])
    return a/len(train_x)
    # return a
def H_1(train_x,train_y,theta):#海森矩阵的逆
    train_x["x3"]=[1.0 for i in range(len(train_x))]
    a = np.zeros((3,3))
    # train_x = pd.concat([DataFrame([1 for i in range(len(train_x))]), train_x], 1)
    for i in range(len(train_x)):
        #a = a-np.array(h(theta, np.mat(train_x.iloc[i, :]).T)*(1-h(theta, np.mat(train_x.iloc[i, :]).T)))[0][0]*np.mat(train_x.iloc[i, :]).T*np.mat(train_x.iloc[i, :])
        a = a + np.eye(3,3)*(np.array(np.exp(np.dot(theta,np.mat(train_x.iloc[i, :]).T)))[0][0]/(np.array(np.exp(np.dot(theta,np.mat(train_x.iloc[i, :]).T)))[0][0]+1)**2)* np.mat(train_x.iloc[i, :]).T*np.mat(train_x.iloc[i, :])
    return (a / len(train_x)).I
    #return a.I
def para(train_x,train_y,theta,times):#输出迭代后的theta,times为限制迭代次数
    for i in range(times):
        itheta=j(train_x,train_y,theta)*H_1(train_x,train_y,theta)
        theta=theta-itheta/abs(min(np.array(itheta)[0]))
    return theta
def pic(train_x,train_y,theta,times):#画图
    theta=para(train_x,train_y,theta,times)
    plot_x=[1,10]
    plot_y=[((-1)*np.array(theta)[0][2]+(-1)*np.array(theta)[0][0]*plot_x[0])/np.array(theta)[0][1],((-1)*np.array(theta)[0][2]+(-1)*np.array(theta)[0][0]*plot_x[1])/np.array(theta)[0][1]]
    plt.scatter(train_x['x1'], train_x['x2'], c=["r" if i==1 else "b" for i in train_y["y"]])
    plt.plot(plot_x, plot_y)

结果如图所示：

 

　　附数据，见网盘链接：http://pan.baidu.com/s/1qXTeDZI 密码：iwuh