从机器学习到逻辑回归
今天，我们只关注机器学习到线性回归这条线上的概念。别的以后再说。为了让大家听懂，我这次也不查维基百科了，直接按照自己的理解用大白话说，可能不是很严谨。
机器学习就是机器可以自己学习，而机器学习的方法就是利用现有的数据和算法，解出算法的参数。从而得到可以用的模型。
监督学习就是利用已有的数据（我们叫X，或者特征），和数据的标注（我们叫Y），找到x和y之间的对应关系，或者说是函数f。
回归分析是一种因变量为连续值得监督学习。而分类是一种应变量为非连续值的监督学习。
这里顺便提一句非连续值和连续值的英文有很多表述。
连续值可以是continuous, numerical, quantitative等。
非连续值可以是categorical, nominal, qualitative等。
逻辑回归虽然名字里面有回归两个字，但是它是分类分析，不是回归分析。逻辑回归，它之所以叫这个名字，是因为它和线性回归实在是太像了。
问题
这里，我们使用sklearn自带的癌症数据集。首先读入数据并放入pandas里面。
import numpy as npimport pandas as pdfrom matplotlib import pyplot as plt%matplotlib inlinefrom sklearn.datasets import load_breast_cancer#from sklearn.linear_model import LogisticRegressionfrom sklearn.model_selection import train_test_splitdataset = load_breast_cancer()data = pd.DataFrame(data=dataset.data, columns=dataset.feature_names)data['cancer'] = [dataset.target_names[t] for t in dataset.target]
输出两个分类。
print(dataset.target_names)
[‘malignant’ ‘benign’]
翻译成中文是[‘恶性’ ‘良性’]。
输出
data[18:28]
这里一共有30个属性。
手写算法
无论是线性回归，逻辑回归，以及我以后会写文章的神经网络，他们的基本思路都是一样的。首先，构建一个函数模型，用这个函数表示从x到y的映射关系。
然后构建一个损失函数loss function。它描述了模型函数f(x)f(x)f(x)和真实值yyy之间的差距。当然，这个差距越小越好。
最后是优化方法。优化方法首先计算损失函数对参数的导数。既，损失函数随参数是变大还是变小的。如果损失函数随着参数的变大而变小，则说明参数应该变大。如果损失函数随着参数的变小而变小，则说明参数应该变小。优化方法里面的α\alphaα是学习率，一个小于1的数值，可以是0.01，0.001， 甚至更小。
模型函数
这里，我们的y值，并非连续的。要么y=0y=0y=0，要么y=1y=1y=1。所以，和线性回归相比，我们要把yyy控制在0和1之间。这时，前人引进了sigmoid函数。当x大于0时，y无限接近于1；当x小于0时，y无限接近于0；当x等于0时，y=0.5。
y^=11+e−z\hat{y}=\frac{1}{1 + e^{- z}}y^=1+e−z1 其中 z=θxz=\theta xz=θx
def sigmoid(z):s = 1/(1+np.exp(-z))s = s.reshape(s.shape[0],1)return s
我们可以把这个函数画出来看看。
def draw_sigmoid():x = np.arange(-6, 6, .01)y = sigmoid(x)plt.plot(x, y, color='red', lw=2)plt.show()draw_sigmoid()
最终，我们的模型函数是：
def model(theta, X):z = np.sum(theta.T * X, axis=1)return sigmoid(z)
损失函数
这里，损失函数用的是cross entropy，的定义如下。
J=−y∗log(y^)−(1−y)∗log(1−y^)J= - y * log(\hat{y}) - (1-y) * log(1-\hat{y})J=−y∗log(y^)−(1−y)∗log(1−y^)
#cross_entropydef cross_entropy(y, y_hat):n_samples = y.shape[0]return sum(-y*np.log(y_hat)-(1-y)*np.log(1-y_hat))/n_samplesdef cost_function(theta, X, y):y_hat = model(theta, X)return cross_entropy(y, y_hat)
优化函数
这里先要解决∂J∂θ\frac{\partial J}{\partial \theta}∂θ∂J
因为有
J=−y∗log(y^)−(1−y)∗log(1−y^)J= - y * log(\hat{y}) - (1-y) * log(1-\hat{y})J=−y∗log(y^)−(1−y)∗log(1−y^)
所以
∂J∂θ=−∂(y∗log(y^)+(1−y)∗log(1−y^))∂θ\frac{\partial J}{\partial \theta} = - \frac{\partial (y * log(\hat{y}) + (1-y) * log(1-\hat{y}))}{\partial \theta}∂θ∂J=−∂θ∂(y∗log(y^)+(1−y)∗log(1−y^))
=−∂(y∗log(y^))∂θ−∂((1−y)∗log(1−y^))∂θ= - \frac{\partial (y * log(\hat{y}))}{\partial \theta} - \frac{\partial ((1-y) * log(1-\hat{y}))}{\partial \theta}=−∂θ∂(y∗log(y^))−∂θ∂((1−y)∗log(1−y^))
=−y∗∂log(y^)∂θ−(1−y)∗log(1−y^)∂θ= - y * \frac{\partial log(\hat{y})}{\partial \theta} - (1-y) * \frac{log(1-\hat{y})}{\partial \theta}=−y∗∂θ∂log(y^)−(1−y)∗∂θlog(1−y^)
=−yy^∗∂y^∂θ−1−y1−y^∗∂(1−y^)∂θ= - \frac{y}{\hat{y}} *\frac{\partial \hat{y}}{\partial \theta} - \frac{1-y}{1-\hat{y}} *\frac{\partial (1-\hat{y})}{\partial \theta}=−y^y∗∂θ∂y^−1−y^1−y∗∂θ∂(1−y^)
=−yy^∗∂y^∂θ+1−y1−y^∗∂y^∂θ= - \frac{y}{\hat{y}} *\frac{\partial \hat{y}}{\partial \theta} + \frac{1-y}{1-\hat{y}} *\frac{\partial \hat{y}}{\partial \theta}=−y^y∗∂θ∂y^+1−y^1−y∗∂θ∂y^
=(−yy^+1−y1−y^)∗∂y^∂θ= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * \frac{\partial \hat{y}}{\partial \theta}=(−y^y+1−y^1−y)∗∂θ∂y^
=(−yy^+1−y1−y^)∗∂(11+e−θx)∂a= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * \frac{\partial (\frac{1}{1 + e^{- \theta x}})} {\partial a}=(−y^y+1−y^1−y)∗∂a∂(1+e−θx1)
=(−yy^+1−y1−y^)∗(−y^2)∗∂(1+e−θx)∂θ= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * (- \hat{y} ^ 2) * \frac{\partial (1 + e^{- \theta x})} {\partial \theta}=(−y^y+1−y^1−y)∗(−y^2)∗∂θ∂(1+e−θx)
=(−yy^+1−y1−y^)∗(−y^2)∗∂(e−θx)∂θ= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * (- \hat{y} ^ 2) * \frac{\partial (e^{- \theta x})} {\partial \theta}=(−y^y+1−y^1−y)∗(−y^2)∗∂θ∂(e−θx)
=(−yy^+1−y1−y^)∗(−y^2)∗e−θx∗∂(−θx)∂θ= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * (- \hat{y} ^ 2) * e^{- \theta x} * \frac{\partial (- \theta x)} {\partial \theta}=(−y^y+1−y^1−y)∗(−y^2)∗e−θx∗∂θ∂(−θx)
=(−yy^+1−y1−y^)∗y^2∗e−θx∗x= (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * \hat{y} ^ 2 * e^{- \theta x} * x=(−y^y+1−y^1−y)∗y^2∗e−θx∗x
∵ y^=11+e−θx\hat{y}=\frac{1}{1 + e^{- \theta x}}y^=1+e−θx1
∴ 1+e−θx=1y^1 + e^{- \theta x}=\frac{1}{\hat{y}}1+e−θx=y^1
∴ e−θx=1y^−1e^{- \theta x}=\frac{1}{\hat{y}} - 1e−θx=y^1−1
∴ e−θx=1−y^y^e^{- \theta x}=\frac{1-\hat{y}}{\hat{y}}e−θx=y^1−y^
∂J∂θ=(−yy^+1−y1−y^)∗y^2∗1−y^y^∗x\frac{\partial J}{\partial \theta} = (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * \hat{y} ^ 2 * \frac{1-\hat{y}}{\hat{y}} * x∂θ∂J=(−y^y+1−y^1−y)∗y^2∗y^1−y^∗x
∂J∂θ=(−yy^+1−y1−y^)∗y^∗(1−y^)∗x\frac{\partial J}{\partial \theta} = (- \frac{y}{\hat{y}} + \frac{1-y}{1-\hat{y}} ) * \hat{y} * (1-\hat{y}) * x∂θ∂J=(−y^y+1−y^1−y)∗y^∗(1−y^)∗x
∂J∂θ=(−y∗(1−y^)+(1−y)∗y^)∗x\frac{\partial J}{\partial \theta} = (- y * (1-\hat{y}) + (1-y) * \hat{y} ) * x∂θ∂J=(−y∗(1−y^)+(1−y)∗y^)∗x
∂J∂θ=(−y+y∗y^+y^−y∗y^)∗x\frac{\partial J}{\partial \theta} = (- y + y *\hat{y} + \hat{y}-y * \hat{y} ) * x∂θ∂J=(−y+y∗y^+y^−y∗y^)∗x
∂J∂θ=(y^−y)∗x\frac{\partial J}{\partial \theta} = (\hat{y}-y ) * x∂θ∂J=(y^−y)∗x
最后，得到优化函数：
θ=θ−α∗∂J∂θ\theta = \theta - \alpha * \frac{\partial J}{\partial \theta}θ=θ−α∗∂θ∂J
θ=θ−α∗(y^−y)∗x\theta = \theta - \alpha * (\hat{y}-y ) * xθ=θ−α∗(y^−y)∗x
python函数如下：
def optimize(theta,X,y):n = X.shape[0]alpha = 1e-1y_hat = model(theta,X)dtheta = (1.0/n) * ((y_hat-y)*X)dtheta = np.sum(dtheta, axis=0)dtheta=dtheta.reshape((31,1))theta = theta - alpha * dthetareturn theta
评估
分类问题的评估比较简单，一般用准确率就可以了。当然也可以用别的（召回率，Precision，ROC，F1 Score等等），其公式如下。
def predict_proba(theta, X):y_hat=model(theta, X)return y_hatdef predict(X, theta):y_hat=predict_proba(theta,X)y_hard=(y_hat > 0.5) * 1return y_harddef accuracy(theta, X, y):y_hard=predict(X, theta)count_right=sum(y_hard == y)return count_right*1.0/len(y)
循环函数
我们不断的调用优化函数来更新参数。
def iterate(theta,X,y,times):costs = []accs = []for i in range(times):theta = optimize(theta,X,y)costs.append(cost_function(theta, X, y))accs.append(accuracy(theta, X, y))return theta, costs, accs
使用算法
准备数据
加载数据
X = dataset.datay = dataset.targetn_features = X.shape[1]
归一化
std=X.std(axis=0)mean=X.mean(axis=0)X_norm = (X - mean) / std
在x矩阵前面加上一列1，这样做的好处是，不需要单独处理截距(interception)。前面的优化方法的推导，已经很麻烦了。如果把截距和斜率(coefficient)分开处理，我又要推导一遍。并且这样程序变得很复杂，很难维护了。
def add_ones(X):ones=np.ones((X.shape[0],1))X_with_ones=np.hstack((ones, X))return X_with_onesX_with_ones = add_ones(X_norm)
求参数
首先，初始化参数。我这里都为1。
theta = np.ones((n_features+1,1))
接着，就可以一直循环求参数了。
theta, costs, accs = iterate(theta, X_train, y_train, 1500)
随着不断的训练，loss不断下降。
（放大到1到100）
而准确率不断提升。
（放大到1到100）
我们最终的loss和准确率是
print(costs[-1], accs[-1])
[0.0489982] [0.99246231]
用测试数据评估为
accuracy(theta, X_test, y_test)
array([0.97660819])
用sklearn计算
from sklearn.linear_model import LogisticRegressionX = dataset.datay = dataset.targetX_train, X_test, y_train, y_test = train_test_split(X_with_ones, y, test_size = 0.3, random_state=12345)lr = LogisticRegression()lr.fit(X_train, y_train)
训练数据集准确率
lr.score(X_train, y_train)
0.992462311557789
测试数据集准确率
lr.score(X_test, y_test)
0.9766081871345029
可以看到，这里的准确率和我们手写的逻辑回归是一样的。大功告成了！
模型参数解释
把函数
y^=11+e−θx\hat{y}=\frac{1}{1 + e^{- \theta x}}y^=1+e−θx1
分解得到
y^=11+e−θ0−θ1x1−θ2x2−θ3x3\hat{y}=\frac{1}{1 + e^{- \theta_0 - \theta_1 x_1 - \theta_2 x_2 - \theta_3 x_3}}y^=1+e−θ0−θ1x1−θ2x2−θ3x31
y^=11+e−θ0e−θ1x1e−θ2x2e−θ3x3\hat{y}=\frac{1}{1 + e^{- \theta_0} e^{- \theta_1 x_1} e^{ - \theta_2 x_2 } e^{- \theta_3 x_3}}y^=1+e−θ0e−θ1x1e−θ2x2e−θ3x31
当θ1\theta_1θ1为0时，e−θ1x1=1e^{- \theta_1 x_1}=1e−θ1x1=1 ，x1x_1x1的变化对结果没有影响。
当θ1\theta_1θ1为大于0时，e−θ1x1e^{- \theta_1 x_1}e−θ1x1随x1x_1x1的增大而变小 ，而分母变小则y^\hat{y}y^变大。既，y^\hat{y}y^随着x1x_1x1的变大而变大。
反之，当θ1\theta_1θ1为小于0时，y^\hat{y}y^随着x1x_1x1的变大而变小。
我们这里以mean radius为例，它的参数为-0.356072，所以，y^\hat{y}y^随着mean radius的变大而变小。而0对应的是恶性，1对应的是良性。我们可以说，mean radius越大，越有可能是恶行肿瘤。
python机器学习手写算法系列
完整源代码：
/juwikuang/machine_learning_step_by_step
欢迎阅读本系列其他文章：
《python机器学习手写算法系列——线性回归》
《python机器学习手写算法系列——逻辑回归》
《python机器学习手写算法系列——决策树》
《python机器学习手写算法系列——kmeans聚类》
