今天，我们只关注机器学习到线性回归这条线上的概念。别的以后再说。为了让大家听懂，我这次也不查维基百科了，直接按照自己的理解用大白话说，可能不是很严谨。
机器学习就是机器可以自己学习，而机器学习的方法就是利用现有的数据和算法，解出算法的参数。从而得到可以用的模型。
监督学习就是利用已有的数据（我们叫X，或者特征），和数据的标注（我们叫Y），找到x和y之间的对应关系，或者说是函数f。
回归分析是一种因变量为连续值得监督学习。而分类是一种应变量为非连续值的监督学习。
这里顺便提一句非连续值和连续值的英文有很多表述。
连续值可以是continuous, numerical, quantitative等。
非连续值可以是categorical, nominal, qualitative等。
逻辑回归虽然名字里面有回归两个字，但是它是分类分析，不是回归分析。逻辑回归，它之所以叫这个名字，是因为它和线性回归实在是太像了。
这里，我们使用sklearn自带的癌症数据集。首先读入数据并放入pandas里面。
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
%matplotlib inline
from sklearn.datasets import load_breast_cancer
#from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
dataset = 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]
| no | mean radius | mean texture | mean perimeter | mean area | mean smoothness | mean compactness | mean concavity | mean concave points | mean symmetry | mean fractal dimension | … | worst texture | worst perimeter | worst area | worst smoothness | worst compactness | worst concavity | worst concave points | worst symmetry | worst fractal dimension | cancer |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 18 | 19.810 | 22.15 | 130.00 | 1260.0 | 0.09831 | 0.10270 | 0.14790 | 0.09498 | 0.1582 | 0.05395 | … | 30.88 | 186.80 | 2398.0 | 0.1512 | 0.3150 | 0.53720 | 0.23880 | 0.2768 | 0.07615 | malignant |
| 19 | 13.540 | 14.36 | 87.46 | 566.3 | 0.09779 | 0.08129 | 0.06664 | 0.04781 | 0.1885 | 0.05766 | … | 19.26 | 99.70 | 711.2 | 0.1440 | 0.1773 | 0.23900 | 0.12880 | 0.2977 | 0.07259 | benign |
| 20 | 13.080 | 15.71 | 85.63 | 520.0 | 0.10750 | 0.12700 | 0.04568 | 0.03110 | 0.1967 | 0.06811 | … | 20.49 | 96.09 | 630.5 | 0.1312 | 0.2776 | 0.18900 | 0.07283 | 0.3184 | 0.08183 | benign |
| 21 | 9.504 | 12.44 | 60.34 | 273.9 | 0.10240 | 0.06492 | 0.02956 | 0.02076 | 0.1815 | 0.06905 | … | 15.66 | 65.13 | 314.9 | 0.1324 | 0.1148 | 0.08867 | 0.06227 | 0.2450 | 0.07773 | benign |
| 22 | 15.340 | 14.26 | 102.50 | 704.4 | 0.10730 | 0.21350 | 0.20770 | 0.09756 | 0.2521 | 0.07032 | … | 19.08 | 125.10 | 980.9 | 0.1390 | 0.5954 | 0.63050 | 0.23930 | 0.4667 | 0.09946 | malignant |
| 23 | 21.160 | 23.04 | 137.20 | 1404.0 | 0.09428 | 0.10220 | 0.10970 | 0.08632 | 0.1769 | 0.05278 | … | 35.59 | 188.00 | 2615.0 | 0.1401 | 0.2600 | 0.31550 | 0.20090 | 0.2822 | 0.07526 | malignant |
| 24 | 16.650 | 21.38 | 110.00 | 904.6 | 0.11210 | 0.14570 | 0.15250 | 0.09170 | 0.1995 | 0.06330 | … | 31.56 | 177.00 | 2215.0 | 0.1805 | 0.3578 | 0.46950 | 0.20950 | 0.3613 | 0.09564 | malignant |
| 25 | 17.140 | 16.40 | 116.00 | 912.7 | 0.11860 | 0.22760 | 0.22290 | 0.14010 | 0.3040 | 0.07413 | … | 21.40 | 152.40 | 1461.0 | 0.1545 | 0.3949 | 0.38530 | 0.25500 | 0.4066 | 0.10590 | malignant |
| 26 | 14.580 | 21.53 | 97.41 | 644.8 | 0.10540 | 0.18680 | 0.14250 | 0.08783 | 0.2252 | 0.06924 | … | 33.21 | 122.40 | 896.9 | 0.1525 | 0.6643 | 0.55390 | 0.27010 | 0.4264 | 0.12750 | malignant |
| 27 | 18.610 | 20.25 | 122.10 | 1094.0 | 0.09440 | 0.10660 | 0.14900 | 0.07731 | 0.1697 | 0.05699 | … | 27.26 | 139.90 | 1403.0 | 0.1338 | 0.2117 | 0.34460 | 0.14900 | 0.2341 | 0.07421 | malignant |
这里一共有30个属性。
无论是线性回归，逻辑回归，以及我以后会写文章的神经网络，他们的基本思路都是一样的。首先，构建一个函数模型，用这个函数表示从x到y的映射关系。
然后构建一个损失函数loss function。它描述了模型函数 f ( x ) f(x) f(x)和真实值 y y y之间的差距。当然，这个差距越小越好。
最后是优化方法。优化方法首先计算损失函数对参数的导数。既，损失函数随参数是变大还是变小的。如果损失函数随着参数的变大而变小，则说明参数应该变大。如果损失函数随着参数的变小而变小，则说明参数应该变小。优化方法里面的 α \alpha α是学习率，一个小于1的数值，可以是0.01，0.001， 甚至更小。
这里，我们的y值，并非连续的。要么 y = 0 y=0 y=0，要么 y = 1 y=1 y=1。所以，和线性回归相比，我们要把 y y y控制在0和1之间。这时，前人引进了sigmoid函数。当x大于0时，y无限接近于1；当x小于0时，y无限接近于0；当x等于0时，y=0.5。
y ^ = 1 1 + e − z \hat{y}=\frac{1}{1 + e^{- z}} y^=1+e−z1 其中 z = θ x z=\theta x z=θ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 ∗ l o g ( y ^ ) − ( 1 − y ) ∗ l o g ( 1 − y ^ ) J= - y * log(\hat{y}) - (1-y) * log(1-\hat{y}) J=−y∗log(y^)−(1−y)∗log(1−y^)
#cross_entropy
def 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_samples
def 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 ∗ l o g ( y ^ ) − ( 1 − y ) ∗ l o g ( 1 − y ^ ) J= - y * log(\hat{y}) - (1-y) * log(1-\hat{y}) J=−y∗log(y^)−(1−y)∗log(1−y^)
所以
∂ J ∂ θ = − ∂ ( y ∗ l o g ( y ^ ) + ( 1 − y ) ∗ l o g ( 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 ∗ l o g ( y ^ ) ) ∂ θ − ∂ ( ( 1 − y ) ∗ l o g ( 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 ∗ ∂ l o g ( y ^ ) ∂ θ − ( 1 − y ) ∗ l o g ( 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^)
= − y y ^ ∗ ∂ y ^ ∂ θ − 1 − y 1 − 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^)
= − y y ^ ∗ ∂ y ^ ∂ θ + 1 − y 1 − 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^
= ( − y y ^ + 1 − y 1 − 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^
= ( − y y ^ + 1 − y 1 − y ^ ) ∗ ∂ ( 1 1 + 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)
= ( − y y ^ + 1 − y 1 − 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)
= ( − y y ^ + 1 − y 1 − 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)
= ( − y y ^ + 1 − y 1 − 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)
= ( − y y ^ + 1 − y 1 − 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 ^ = 1 1 + e − θ x \hat{y}=\frac{1}{1 + e^{- \theta x}} y^=1+e−θx1
∴ 1 + e − θ x = 1 y ^ 1 + e^{- \theta x}=\frac{1}{\hat{y}} 1+e−θx=y^1
∴ e − θ x = 1 y ^ − 1 e^{- \theta x}=\frac{1}{\hat{y}} - 1 e−θx=y^1−1
∴ e − θ x = 1 − y ^ y ^ e^{- \theta x}=\frac{1-\hat{y}}{\hat{y}} e−θx=y^1−y^
∂ J ∂ θ = ( − y y ^ + 1 − y 1 − 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 ∂ θ = ( − y y ^ + 1 − y 1 − 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-1
y_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 * dtheta
return theta
分类问题的评估比较简单，一般用准确率就可以了。当然也可以用别的（召回率，Precision，ROC，F1 Score等等），其公式如下。
def predict_proba(theta, X):
y_hat=model(theta, X)
return y_hat
def predict(X, theta):
y_hat=predict_proba(theta,X)
y_hard=(y_hat > 0.5) * 1
return y_hard
def 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.data
y = dataset.target
n_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_ones
X_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])
from sklearn.linear_model import LogisticRegression
X = dataset.data
y = dataset.target
X_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 ^ = 1 1 + e − θ x \hat{y}=\frac{1}{1 + e^{- \theta x}} y^=1+e−θx1
分解得到
y ^ = 1 1 + e − θ 0 − θ 1 x 1 − θ 2 x 2 − θ 3 x 3 \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 ^ = 1 1 + e − θ 0 e − θ 1 x 1 e − θ 2 x 2 e − θ 3 x 3 \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 − θ 1 x 1 = 1 e^{- \theta_1 x_1}=1 e−θ1x1=1 ， x 1 x_1 x1的变化对结果没有影响。
当 θ 1 \theta_1 θ1为大于0时， e − θ 1 x 1 e^{- \theta_1 x_1} e−θ1x1随 x 1 x_1 x1的增大而变小 ，而分母变小则 y ^ \hat{y} y^变大。既， y ^ \hat{y} y^随着 x 1 x_1 x1的变大而变大。
反之，当 θ 1 \theta_1 θ1为小于0时， y ^ \hat{y} y^随着 x 1 x_1 x1的变大而变小。
我们这里以mean radius为例，它的参数为-0.356072，所以， y ^ \hat{y} y^随着mean radius的变大而变小。而0对应的是恶性，1对应的是良性。我们可以说，mean radius越大，越有可能是恶行肿瘤。
完整源代码：
https://github.com/juwikuang/machine_learning_step_by_step
欢迎阅读本系列其他文章：