对率回归的实验

对数几率回归在python中的实现

在做分类任务时，需要找一个单调可微函数将分类任务的真实标记y与线性回归模型所预测的值联系起来。对数几率函数是用来“替代”单位跃阶函数的，满足单调可微的条件。以下是对数几率函数：$y=1/{(1+e^{-z})}$，其中$z= w^Tx+b$。对其进行取对数，得到

$$ln(y/(1-y))=w^Tx+b$$
该公式的解释：将y视为样本x作为正样本的可能性，1-y视为其反列的可能性，对几率取对数得到对数几率。
若将y视为类后验概率估计p(y=1| x)，则对数几率可以写为： $$ln(p(y=1|x)/p(y=0|x))=w^Tx+b$$
通过极大似然法来估计 w和b。对给定的数据{{ $x_i$, $y_i$}} $^m_{i=1}$，对数回归模型最大化‘对数似然 $$l(\mathbf{w},b)=\sum_{i=1}^m lnp(y_i|\mathbf{x_i;w,}b)$$
改写似然项后可以得到最终公式为 $$l(\beta) = \sum_{i=1}^m (-y_i\beta^T\mathbf x_i + ln(1+exp(\mathbf{\beta^Tx_i})))$$
求得最优解 $\beta ^* = arg min_\beta l(\beta)$
本文中采用牛顿法：
$$\beta ^{t+1} = \beta ^t - (\frac{\partial^{2}l(\beta)}{\partial \beta \partial \beta^T})^-{1} (\frac{\partial l(\beta)}{\partial \beta})$$
分别求得一阶导数和二阶导数得到：
$$\frac{\partial l(\beta)}{\partial \beta} = -\sum_{i=1}^m\mathbf x_i(y_-p_1(\mathbf {x_i;\beta}))$$
$$\frac{\partial^{2}l(\beta)}{\partial \beta \partial \beta^T} = \sum_{i=1}^m \mathbf {x_i x_i^T}p_1(\mathbf{ x_i;\beta})(1-p_1(\mathbf{x_i;\beta}))$$

以上就是对率回归的模型。以下在《机器学习》书中的西瓜数据3.0中实验：

# -*- coding: utf-8 -*-
import numpy as np
import pandas as pd

inputfile = 'xigua3.0.xlsx'

#导入数据
data = pd.read_excel(inputfile, 'Sheet1')
x = np.array([list(data[u'密度']),list(data[u'含糖率']),[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]])
x = x.T
y = np.array([1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0])

#初始化参数
beta = np.array([[0.01],[0.03],[1]])
l_beta = 0
old_l_beta = 0
n = 0
while True:
    #计算l(beta)

    beta_T = np.transpose(beta)
   # print(np.array([x[0,:]]).shape)
    for i in np.arange(len(x)):
        l_beta = l_beta + (-y[i]*np.dot(beta_T, np.array([x[i,:]]).T) +  np.log(1+np.exp(np.dot(beta_T,np.array([x[i,:]]).T))))

    if np.abs(l_beta - old_l_beta).any()<=0.000001:
        break
    #进行迭代
    dbeta = 0
    d2beta = 0
    n = n+1
    old_l_beta = l_beta
    for i in np.arange(len(x)):
        x_i = np.array([x[i,:]])
        x_i_2 = np.dot(x_i,x_i.T)
        exp_b_x = np.exp(np.dot(np.transpose(beta),x_i.T))  

        dbeta = dbeta - np.array([x[i,:]])*( y[i]-( exp_b_x/(1+exp_b_x)))

        d2beta = d2beta + x_i_2*exp_b_x/((1+exp_b_x)*(1+exp_b_x))
    beta = beta - np.dot(np.linalg.inv(d2beta),dbeta).T

    print("迭代次数=>", n)
    print('模型参数=>', beta)
if __name__ == '__main__':
    print(beta)