python 逻辑回归demo

Logistic回归–未正则化

回归中用到的公式列举如下:

sigmoid 函数

用线性回归解决分类问题是不妥的，因此需要对其做一个转化，令：

$h_\theta (x)=g(\theta^TX)$

g 代表一个常用的逻辑函数（logistic function）为S形函数（Sigmoid function），公式为：$g\left( z \right)=\frac{1}{1+{{e}^{-z}}}$

合起来，我们得到逻辑回归模型的假设函数：

${{h}_{\theta }}\left( x \right)=\frac{1}{1+{{e}^{-{{\theta }^{T}}X}}}$

本质上将一个线性表达式映射到了sigmoid函数中，该函数是一个S型曲线，当 $\theta^TX>0$ 时，$g(z)\geq0.5$ ，将其分类为$1$，否则输出$0$。具体如下图：

cost function(损失函数)

• $max(\ell(\theta)) = min(-\ell(\theta))$
• choose $-\ell(\theta)$ as the cost function
  这里没有用线性回归的损失函数是因为只有当损失函数为凸函数的时候才能通过求导得出极值点从而输出局部最优解，而逻辑回归里的假设函数$h_\theta (x)$是一个非线性函数，不适用于线性回归的代价函数，因此我们对其取对数化，将其构造成一个$Cost(h_\theta (x),y)$取值为$[0 , \infty)$的函数，函数具体表现如下：

$J\left( \theta \right)=-\frac{1}{m}\sum\limits_{i=1}^{m}{[{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)+\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}$
$=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}$

gradient descent(梯度下降)

• 批量梯度下降（batch gradient descent）
• 向量化计算： $\frac{1}{m} X^T( Sigmoid(X\theta^T) - y )$

$\frac{\partial J\left( \theta \right)}{\partial {{\theta }_{j}}}=\frac{1}{m}\sum\limits_{i=1}^{m}{({{h}_{\theta }}\left( {{x}^{(i)}} \right)-{{y}^{(i)}})x_{_{j}}^{(i)}}$

下面开始构建该回归模型

'''准备包包'''
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('fivethirtyeight') #样式美化
from sklearn.metrics import classification_report#这个包是评价报告

'''准备数据'''
data = pd.read_csv('ex2data1.txt', names=['exam1', 'exam2', 'admitted'])
data.head()
data.describe()
sns.set(context="notebook", style="darkgrid", palette=sns.color_palette("RdBu", 2),color_codes = False)

'''瞅瞅长啥样'''
sns.lmplot('exam1', 'exam2', hue='admitted', data=data, 
           size=6, 
           fit_reg=False, #控制是否显示拟合直线
           scatter_kws={"s": 50} #s代表圆形，50代表点的大小
          )
plt.show()#看下数据的样子

def get_X(df):#读取特征
	'''
	这里用的是concat，和之前的insert是一样的作用
	'''
    #data = df.insert(0,'ones',1) #和下面的concat函数可以实现同样的功能
    ones = pd.DataFrame({'ones': np.ones(len(df))})#ones是m行1列的dataframe
    data = pd.concat([ones, df], axis=1)  # 合并数据，根据列合并
    return data.iloc[:, :-1].as_matrix()  # 这个操作返回 ndarray,不是矩阵，之后的版本不再用as_matrix()，而用.values替代


def get_y(df):#读取标签
 	'''假定最后一列是分类标签'''
    return np.array(df.iloc[:, -1])#df.iloc[:, -1]是指df的最后一列


def normalize_feature(df): #特征归一化
   	'''
    多个特征变量需要做归一化处理，利用apply函数，默认axis=0，即按纵向进行运算，在这里就是每列的均值和标准差，关于apply函数的使用感想，单独查阅
    其中apply函数是对dataframe类型数据进行按列或者按行计算的一种函数，和lambda函数一起用   
    '''
    return df.apply(lambda x: (x - x.mean()) / x.std())#特征缩放

'''将特征值和标签选出'''
X = get_X(data)
print(X.shape)

y = get_y(data)
print(y.shape)

'''定义一个sigmoid函数'''
def sigmoid(z):
    # your code here  (appro ~ 1 lines)
    gz = 1/(1+np.exp(-z))
    return gz    
#如果z是一个矩阵的话，矩阵里面的每个元素也会经历sigmoid函数的变化

'''定义损失函数'''
theta=np.zeros(X.shape[1]) # X(m*n) so theta is n*1
#theta =np.mat(theta)
def cost(theta, X, y):
    ''' cost fn is -l(theta) for you to minimize'''
    costf = np.mean(-y*np.log(sigmoid(X@theta))-(1-y)*np.log(1-sigmoid(X@theta)))
    return costf
cost(theta, X, y)
# X @ theta与X.dot(theta)等价
#dot和*以及@的乘法区别。当两个数据都是数组时，*代表元素相乘，@和.dot为矩阵乘法，后期补充

'''批量梯度下降'''
def gradient(theta, X, y):
    grad = X.T@(sigmoid(X@theta)-y)/len(X)
    return grad
gradient(theta, X, y)

'''拟合参数'''
import scipy.optimize as opt
res = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient) #fun表是选取寻找最低值的目标函数，args是常数值，jac为雅各比行列式方法，选取自己写的梯度下降
print(res)

fun: 0.20349770159041028
jac: array([-1.36170577e-07, -2.02377035e-05, 6.14692479e-06])
message: ‘Optimization terminated successfully.’
nfev: 72
nhev: 0
nit: 28
njev: 243
status: 0
success: True
x: array([-25.16141434, 0.20623236, 0.20147225])

'''将预测变量输出为0，1'''
def predict(x, theta):
    y_pred = sigmoid(X@theta)
    return (y_pred>=0.5).astype(int)#布尔值转换
'''预测分类效果验证'''
final_theta = res.x
y_pred = predict(X, final_theta)
print(classification_report(y, y_pred))

输出结果如下：

precision    recall  f1-score   support

       0       0.87      0.85      0.86        40
       1       0.90      0.92      0.91        60

accuracy                           0.89       100

'''在图形上表示画出决策边界，即找到theta^TX=0的点'''
coef = -(res.x / res.x[2])  # res.x是一个一行三列的数组，分别代表theta0,theta1,theta2，这里全部除以theta2，其实除以theta1也行，将其表现成y=-theta0/theta2 - theta1/theta2 *x 的形式
print(coef)
x = np.arange(130, step=0.1)
y = coef[0] + coef[1]*x

'''画图'''
sns.set(context="notebook", style="ticks", font_scale=1.5)
sns.lmplot('exam1', 'exam2', hue='admitted', data=data, 
           size=6, 
           fit_reg=False, 
           scatter_kws={"s": 25}
          )

plt.plot(x, y, 'grey')
plt.xlim(0, 130)
plt.ylim(0, 130)
plt.title('Decision Boundary')
plt.show()