机器学习

补遗: 机器学习手记系列 3 线性回归样例程序

距离 2012 的两三年后(这篇的草稿时间)又过了两三年,这个补遗看起来也烂尾了 -.-

之前在机器学习手记系列 3: 线性回归和最小二乘法后面留了个问题, 也给了结果, 但是当时说好的程序代码并没给出来, 那个手记系列的坑感觉填不上了, 但是已经刨好的小坑还是填上吧

现在已经有很多深度学习框架和教程来教这个,自己也忘得差不多了,就不班门弄斧裸写。推荐看一下 动手学深度学习 http://zh.gluon.ai/index.html,Deep Learning 领域大神 李沐 等人在维护(我能凑不要脸的蹭热度说下这是前百度同事我们还一起吃饭打牌来着么)。刨的小坑就按 线性回归的从零开始实现 http://zh.gluon.ai/chapter_deep-learning-basics/linear-regression-scratch.html 里的做法来实现

先重复下问题

如下式子里不同的阿拉伯数字只是一个符号, 实际表示的可能是其他数字
967621 = 3
797321 = 1
378581 = 4
422151 = 0
535951 = 1
335771 = 0

根据上述式子, 判断下式等于?
565441 = ?

这题的脑筋急转弯版本答案是看每个数字有几个圈,就代表几,这样 1/2/3/4/5/7 都是 0 个圈,6/9 是 1 个圈,8 是 2 个圈,所以最后 565441 里面只有 6 有 1 个圈,答案为 1

按 gluon 上的教程我们也来走一遍,装环境什么的就看 gluon 了,先引入要用的包

from mxnet import autograd, nd

真正做线性回归是没法只用这么一点数据来模拟的,所以我们要先根据真实值来构造一些数据(这里跟 gluon 不一样的是我没有 bias 因子 b,后面也请一并注意)

num_inputs = 9          # 特征数,当前问题里的变量数 1-9
num_examples = 1000     # 样例数,我们会随机生成多少份样例来学习
true_w = nd.array([0, 0, 0, 0, 0, 1, 0, 2, 1])  # 真实值
features = nd.random.normal(scale=1, shape=(num_examples, num_inputs))  # 随机生成数据集
labels = nd.dot(features, true_w)                                   # 数据集对应的结果

初始化模型参数并创建梯度

w = nd.random.normal(scale=0.01, shape=(9, 1))
w.attach_grad()

定义模型,我们就是做的矩阵乘法

def linreg(X, w):
    return nd.dot(X, w)

定义损失函数,用平方损失

def squared_loss(y_hat, y):
    return (y_hat - y.reshape(y_hat.shape)) ** 2 / 2

定义优化算法,用小批量随机梯度下降(因为我们只用了一个大参数 w,所以还是比 gluon 的样例简单)

def sgd(param, lr, batch_size):
    param[:] = param - lr * param.grad / batch_size

训练,取步长 lr 为 0.01,轮次为 1000 轮

def train():
    lr = 0.01
    num_epochs = 1000
    net = linreg
    loss = squared_loss

    for epoch in range(num_epochs):
        with autograd.record():
            l = loss(net(features, w), labels)
        l.backward()
        sgd(w, lr, labels.size)
        train_l = loss(net(features, w), labels)
        if epoch % 100 == 99:
            print("epoch {}, loss {}, w {}".format(epoch + 1, train_l.mean().asnumpy(), w))

验证下结果看看

if __name__ == "__main__":
    train()
    test = nd.array([1, 0, 0, 2, 2, 1, 0, 0, 0])    # 测试集,565441
    print(nd.dot(test, w))

随便跑了一次输出如下,注意模型里每个值的科学计数法的指数

epoch 1000, loss [  5.72006487e-09], w
[[ -6.20802666e-06]
 [  1.62000088e-05]
 [ -1.03610901e-05]
 [  7.82768348e-06]
 [  2.59973749e-05]
 [  9.99964714e-01]
 [  1.86312645e-05]
 [  1.99990368e+00]
 [  1.00001490e+00]]
<NDArray 9x1 @cpu(0)>

[ 1.00002611]
<NDArray 1 @cpu(0)>

忽略精度问题,可以认为符合真实结果

全部代码详见 https://gist.github.com/whusnoopy/af0aa6fd276ace8a7c4d483e586e936d

机器学习手记系列 3: 线性回归和最小二乘法

好几个月没再继续, 挖坑不填是不对的, 还是继续写手记.

线性回归

线性回归一般用来学习一维自变量和一维因变量之间的线性关系. 如果存在一维自变量 x, 同时还有一维因变量 y = f(x), 如果有一堆对不同 x 下 y 的观测值, 即 的观测对, 且如果 x, y 之间存在较明显的线性关系, 可以用 f(x) = a*x + b 这样的方程表示, 则可以用线性回归的方法学习出 a 和 b 的值, 同时估计这个拟合方法的误差 r.

扩展一下, 线性回归也可以指 y = a1*x1 + a2*x2 + ... + ak*xk + b 这样多维自变量和一维因变量之间的线性关系 (多项式里自变量的最高幂次都是 1), 同样也可以用回归的方式学出来里面不同的系数和常数项的值. 只有一维自变量的称之为一元线性回归, 否则是多元线性回归.

判断线性回归好坏, 一般就用平方误差和来描述, 其表达为 (f(x1)-y1)2 + (f(x2)-y2)2 + ... + (f(xk)-yk)2), 此值如果为 0 则说明自变量和因变量存在完全的线性关系, 否则是近似线性, 越小近似的越好. 这个东西看着有没有一点面熟? 其实就是机器学习手记系列 2: 离线效果评估里最后提到的 MSE 的非平均版本.

解方程法

假如输入样本存在绝对的线性关系, 即最后的误差为 0, 则问题变为解二元一次方程 (一元线性回归里的系数加常数) 或 N+1 元一次方程 (多元线性回归里 N 个自变量的系数加常数). 这没什么好说的, 直接对原输入直接求解就行了, 类似计算方法这样的课本上有的是解法, 做梯度下降或牛顿法乃至矩阵分解都可以解.

梯度下降法

考虑到绝大部分情况下不存在绝对的线性关系, 则问题可以变成怎么求平方误差和的最小值点. 如果是一元线性回归, 目标函数变为 g(a, b) = (f(x1)-y1)2 + (f(x2)-y2)2 + ... + (f(xk)-yk)2

我们的目标就是让这个目标函数值最小, 选定一组 的初始值, 然后求其梯度方向, 每次前进一个小步长, 再求梯度前进, 直到目标函数值不再下降, 说明我们已经走到了一个极值点附近, 终止迭代. 对一元线性回归, 梯度方向是对函数求偏导得到的向量方向.

另外需要注意的是, 梯度下降不一定能找到最优解, 可能会在某个局部最优解那陷进去就出不来了.

这部分更详细的推导请见参考资料里 “机器学习中的数学(1)” 一篇, 里面的公式和图做的很赞, 思路也比我清晰.

牛顿法

梯度下降一般遇到的问题迭代步长不好选, 选太小到极值点太慢, 搞太大又会在极值点附近时因为步长太大跳过去了.

牛顿法最大的贡献就是同时给出了梯度方向和迭代步长, 几乎是一步到位的求解. 方法同解方程一样, 对新的损失目标函数求解, 只是一次解可能还不够好, 需要多做几次迭代. 一般梯度下降可能需要上千轮的迭代, 而牛顿法几次迭代就能到极值点了.

最小二乘法

伟大的高斯同学提出并证明了最小二乘法是最好解答, 证明过程略… 直接看维基或百度百科上的原文吧 (数学不好伤不起).

应用

虽然这个方法看起来很简单粗暴, 但是很多时候变化确实就是线性的. 比如在很多论文和工业实践中, 大家认为同等质量的广告或搜索结果, 放在从上到下不同的位置上, 其点击率和位置的关系符合线性关系, 即 ctr(rank) = a*rank + b.

在六月的一次随笔杂记里, 提到了这样的问题:

如下式子里不同的阿拉伯数字只是一个符号, 实际表示的可能是其他数字
967621 = 3
797321 = 1
378581 = 4
422151 = 0
535951 = 1
335771 = 0

根据上述式子, 判断下式等于?
565441 = ?

假设每个式子最后做的都是加法, 并把字符 0~9 映射到 x0~x9, 则统计不同字母出现的次数就可以列线性返程, 可以将第一个式子表示为

x0*0 + x1*1 + x2*1 + x3*0 + x4*0 + x5*0 + x6*2 + x7*1 + x8*0 + x9*1 = 3

其他类推. 对这堆式子求解就可以得到不同数字对应的真实数值, 可以得到 565441 = 1. // 具体代码和方法下次给出

参考资料

* Wiki Least squares: http://en.wikipedia.org/wiki/Least_squares
* Wiki Mean Squared Error: http://en.wikipedia.org/wiki/Mean_squared_error
* 中文维基 最小二乘法: http://zh.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E4%BA%8C%E4%B9%98%E6%B3%95
* 百度百科 线性回归: http://baike.baidu.com/view/449540.htm
* 百度百科 最小二乘法: http://baike.baidu.com/view/139822.htm
* 人人上的日志 幼儿园的题目和机器学习的关系: http://blog.renren.com/share/30314/13432269197
* 机器学习中的数学(1): http://www.cnblogs.com/LeftNotEasy/archive/2010/12/05/mathmatic_in_machine_learning_1_regression_and_gradient_descent.html

// 本文开写于 9.24, 拖到 10.26 才马马虎虎完成, 最后那个题的解也没写, 各种错乱后续再修改或补充吧

机器学习手记系列 2: 离线效果评估

上一次说到选特征的一个简单方法, 但是如果真的要评估一个方法或者一类特征的效果, 简单的相似度计算是不够的, 在上线实验之前, 还是需要有一些别的方式来做验证.

我遇到过的大部分机器学习问题, 最终都转成了二分类问题 (概率问题). 最直白的比如 A 是否属于集合 S (某照片中的人脸是否是人物 Z), 排序问题也可以转换为二分类问题, 比如广告点击率或推荐的相关度, 把候选集分为点击/不点击或接受推荐/不接受推荐的二分类概率. 那在上线之前, 可以用过一些分类器性能评估的方法来做离线评估.

分类器的正确率和召回率

前几天在无觅上看到有人分享了一篇 数据不平衡时分类器性能评价之ROC曲线分析, 把这个问题已经讲差不多了, 我这复述一下.

先说混淆矩阵 (confusion matrix). 混淆矩阵是评估分类器可信度的一个基本工具, 设实际的所有正样本为 P (real-Positive), 负样本为 N (real-Negative), 分类器分到的正样本标为 pre-Positive’, 负样本标为 pre-Negetive’, 则可以用下面的混淆矩阵表示所有情况:

              | real-positive       | real-negative
pre-positive' | TP (true positive)  | FP (false positive)
pre-negative' | FN (false negative) | TN (true negative)

通过这个矩阵, 可以得到很多评估指标:

FP rate = FP / N
TP rate = TP / P
Accuracy = (TP + TN) / (P + N)    # 一般称之为准确性或正确性
Precision = TP / (TP + FP)        # 另一些领域的准确性或正确性, 经常需要看上下文来判断
Recall = TP / P                   # 一般称之为召回率
F-score = Precision * Recall

在我接触过的大部分工作中, 大家都在关注 Precision 和 Recall. 同引用原文中提到的, 这样的分类评估性能只在数据比较平衡时比较好用 (正负例比例接近), 在很多特定情况下正负例是明显有偏的 (比如万分之几点击率的显示广告), 那就只能作为一定的参考指标.

分类器的排序能力评估

很多情况下我们除了希望分类器按某个阈值将正负样本完全分开, 同时还想知道候选集中不同条目的序关系. 比如广告和推荐, 首先需要一个基础阈值来保证召回的内容都满足基本相关度, 比如我一大老爷们去搜笔记本维修代理你给我出一少女睫毛膏的广告或推荐关注, 我绝对飙一句你大爷的然后开 AdBlock 屏蔽之. 在保证了基础相关性 (即分类器的正负例分开) 后, 则需要比较同样是正例的集合里, 哪些更正点 (其实说白了就是怎样才收益最大化). 一般来说, 如果分类器的输出是一个正例概率, 则直接按这个概率来排序就行了. 如果最终收益还要通过评估函数转换, 比如广告的 eCPM = CTR*Price, 或推荐里 rev = f(CTR), (f(x) 是一个不同条目的获益权重函数), 那么为了评估序是否好, 一般会再引入 ROC 曲线和 AUC 面积两个指标.

ROC 曲线全称是 Receiver Operating Characteristic (ROC curve), 详细的解释可以见维基百科上的英文词条 Receiver_operating_characteristic 或中文词条 ROC曲线. 我对 ROC 曲线的理解是, 对某个样本集, 当前分类器对其分类结果的 FPR 在 x 时, TPR 能到 y. 如果分类器完全准确, 则在 x = 0 时 y 就能到 1, 如果分类器完全不靠谱, 则在 x = 1 时 y 还是为 0, 如果 x = y, 那说明这个分类器在随机分类. 因为两个都是 Rate, 是 [0, 1] 之间的取值, 所以按此方法描的点都在一个 (0, 0), (1, 1) 的矩形内, 拉一条直线从 (0, 0) 到 (1, 1), 如果描点在这条直线上, 说明分类器对当前样本就是随机分的 (做分类最悲催的事), 如果描点在左上方, 说明当前分类器对此样本分类效果好过随机, 如果在右下方, 那说明分类器在做比随机还坑爹的反向分类. 引用维基百科上的一个图来说明:

ROC 曲线示例

其中 C’ 好于 A (都是正向的), B 是随机, C 是一个反效果 (跟 C’ 沿红线轴对称, 就是说把 C 的结果反过来就得到 C’).

如果我们有足够多的样本, 则对一个分类器可以在 ROC 曲线图上画出若干个点, 把这些点和 (0, 0), (1, 1) 连起来求凸包, 就得到了 AUC 面积 (Area Under Curve, 曲线下面积). 非常明显, 这个凸包的最小下面积是 0.5 (从 (0, 0) 到 (1, 1) 的这条线), 最大是 1.0 (整个矩形面积), AUC 值越大, 说明分类效果越好.

用 ROC 曲线定义的方式来描点计算面积会很麻烦, 不过还好前人给了我们一个近似公式, 我找到的最原始出处是 Hand, Till 在 Machine Learning 2001 上的一篇文章给出 [文章链接]. 中间的推导过程比较繁琐, 直接说我对这个计算方法的理解: 将所有样本按预估概率从小到大排序, 然后从 (0, 0) 点开始描点, 每个新的点是在前一个点的基础上, 横坐标加上当前样本的正例在总正例数中的占比, 纵坐标加上当前样本的负例在总负例数中的占比, 最终的终点一定是 (1, 1), 对这个曲线求面积, 即得到 AUC. 其物理意义也非常直观, 如果我们把负例都排在正例前面, 则曲线一定是先往上再往右, 得到的面积大于 0.5, 说明分类器效果比随机好, 最极端的情况就是所有负例都在正例前, 则曲线就是 (0, 0) -> (0, 1) -> (1, 1) 这样的形状, 面积为 1.0.

同样给一份 C 代码实现:

struct SampleNode {
  double predict_value;
  unsigned int pos_num;
  unsigned int neg_num;
};

int cmp(const void *a, const void *b)
{
   SampleNode *aa = (SampleNode *)a;
   SampleNode *bb = (SampleNode *)b;
   return(((aa->predict_value)-(bb->predict_value)>0)?1:-1);
}

double calcAuc(SampleNode samples[], int sample_num) {
  qsort(samples, sample_num, sizeof(SampleNode), cmp);

  // init all counters
  double sum_pos = 0;
  double sum_neg = 0;
  double new_neg = 0;
  double rp = 0;
  for (int i = 0; i < sample_num; ++i) {
    if (samples[i].neg_num >= 0) {
      new_neg += samples[i].neg_num;
    }

    if (samples[i].pos_num >= 0) {
      // calc as trapezium, not rectangle
      rp += samples[i].pos_num * (sum_neg + new_neg)/2;
      sum_pos += samples[i].pos_num;
    }
    sum_neg = new_neg;
  }

  return rp/(sum_pos*sum_neg);
}

分类器的一致性

如果分类器的概率结果就是最终结果序, 那 AUC 值基本可以当作最终效果来用. 但是实际应用中分类器的结果都要再做函数转换才是最终序, 则评估的时候需要将转换函数也带上去做 AUC 评估才行. 某些应用中这个转换函数是不确定的, 比如广告的价格随时会变, 推荐条目的重要性或收益可能也是另一个计算模型的结果. 在这种情况下, 如果我们可以保证分类器概率和实际概率一致, 让后续的转换函数拿到一个正确的输入, 那么在实际应用中才能达到最优性能.

为了评估分类器概率和实际概率的一致性, 引入 MAE (Mean Absolute Error, 平均绝对误差) 这个指标, 维基百科对应的词条是 Mean_absolute_error. 最终的计算方法很简单, 对样本 i, fi 是预估概率, yi 是实际概率, 则 i 上绝对误差是 ei, 累加求平均就是 MAE:

MAE 公式

MAE 的值域是 [0, +∞), 值越小说明分类器输出和实际值的一致性越好. 我个人认为如果 MAE 和实际概率一样大, 那这个分类器的波动效果也大到让预估近似随机了.

MAE 看起来和标准差有点像, 类似标准差和方差的关系, MAE 也有一个对应的 MSE (Mean Squared Error, 均方差?), 这个指标更多考虑的是极坏情况的影响, 计算比较麻烦, 一般用的也不多, 有兴趣的可以看维基百科上的词条 Mean_squared_error.

MAE 计算太简单, MSE 计算太纠结, 所以都不在这给出代码实现.

机器学习手记系列 1: Pearson 相关系数

系列说明

按合总指示, 给人人的机器学习小组写点科普性质的东西. 其实自己好像一直都没去系统的学过这些东西, 都是野路子乱搞, 这里把过去学的一点东西写出来, 记录一下, 班门弄斧, 欢迎拍砖.

自己接触到的机器学习, 几乎都是在用历史预估未来某事件发生的概率 (广告点击率, 推荐接受度, 等等).

将这个过程细化一下, 首先都是对历史样本提取特征, 将样本转换成用特征序列来描述, 将同类事件合并, 然后通过某种拟合方式去让特征带上合适的权重, 用于描述事件发生的概率, 最后对还未发生的同类事件, 同样将其转换成特征序列, 用学出来的权重转换成预估概率.

这里有两个关键问题, 一是特征选取, 二是拟合还原方法. 特征选取是为了将样本做合理拆分合并, 同质的才划分到一起, 不然就最后的预估还是随机猜. 拟合还原方法是保证在对数据做了合理拆分后, 能将特征的权重拟合到原数据上且能在预估时还原成概率.

问题

对怎么选特征, 似乎从来都没有好的普适性方法, 但是怎么验证选的特征靠不靠谱, 方法倒是挺多. 先抛开选特征的指导方向不说 (说也说不清), 如果我们选出了一类描述特征, 怎么验证其效果?

特征效果验证

最直接粗暴也是终极方案就是直接拿到线上去应用, 好就是好, 不好就是不好. 这是一句彻底的废话, 也是真理… 不过实际操作中很难真的去这么做, 一是如果要完成整个流程会比较耗时和麻烦, 二是没有那么多线上资源拿来实验, 三是如果实验不好带来的负面影响会非常大, 广告会损失收入, 推荐会严重影响用户体验. 所以如果线下没验证的心里有谱, 没人敢直接拍上去实验的, 老板也不会让乱来, 都是钱和能转换成钱的用户啊.

退回到离线验证上, 终极离线验证也还是拿机器学习的产出 (分类树, LR 模型, 或别的什么) 去评估一部分实验数据, 然后看对实验数据的预估结果是否和实验数据的实际表现一致. 这个还是存在耗时耗资源的问题, 后面再说, 先说简单的.

不管是分类树, 还是 LR 或别的 boosting 什么的, 都是希望能找到有区分度的特征, 能将未来不同的数据尽可能划开. 如果特征是一个 0/1 特征, 那数据就应该能明显被分成不一样的两份. 比如在豆瓣, “历史上关注过计算机类别书目的人” 是一个 0/1 特征, 如果拿这个特征来分拆人群, 并评估 “未来是否关注计算机类别书目”, 评估指标的 1 绝大部分都落在区分特征的 1 中, 那说明这是一个非常正相关的区分 (曾经干过某事的人会继续干另一件事), 效果很好, 反之拿去评估 “未来是否关注女性言情小说类别书目”, 很可能评估指标的 1 绝大部分都落在区分特征的 0 里, 那是非常负相关的区分 (曾经干过某事的人不会再干另一件特定的事), 效果也挺好, 但是如果是评估 “未来是否会关注武侠”, 评估指标的 1 是比较均匀的散布在区分特征的 0/1 里, 那就说明这个特征对该评估指标没有区分度, 还不如没有. (这些例子是随便拍脑袋写的, 不保证其正确性)

如果是一维连续特征, 则最后总特征总可以转换成一个一维向量 (连续值可以离散化成整数区间), 跟 0/1 特征一样, 比较这个特征的自变量取值向量和对应到评估数据上的取值向量的相关度就能判断效果好坏 (正相关或负相关都是好的). 一般最简单的是使用 Pearson (皮尔生) 乘积矩相关系数 r 来做是否线性相关的判断, 英文 wiki 上的条目 Pearson product-moment correlation coefficient 对该系数的含义和计算方法有比较详细的说明, 中文翻译比较杂, 百度百科上的是皮尔森相关系数, 微软的翻译是皮尔生相似度. 简要的说, pearson 相关系数是一个 [-1, 1] 之间的实数, 取值越接近 -1 表示特征值和评估值越负线性相关, 越接近 1 表示越正相关, 越接近 0 表示越不相关 (只是线性, 可能会有其他相关的关系).

还是拿豆瓣的例子说, 比如 “历史上关注过计算机类别书的数量” 作为一个人群划分特征, 那这维特征的自变量向量会是 <0, 1, 2, ...>, 为了取值方便, 将其截断到超过十本的也等于 10, 则向量变为 <0, 1, ..., 10>, 如果去评估 “继续关注计算机类别书的概率”, 这个概率取值可能是 <0.0, 0.1, ..., 1.0>, 则其 Pearson 相似度会是 1. (当然, 这个例子举的太假了点, 先这么着吧)

将问题化简为算特征取值和评估指标的 Pearson 相似度后需要做的工作就会少很多, 直接跑个 Map/Reduce 从 LOG 里提取下数据, 然后看看值的相关性就知道特征是否有区分度了, 没区分度的可以先不考虑 (或者将曲线相关的可能转变成线性关系, 再判断), 有区分度的才继续走更完整的离线验证, 加入分类树或概率模型和其他特征一起作用看效果, 如果还不错就上线实验.

微软的 Excel 里就带了 Pearson 相似度的计算公式, 可以很方便的拿来评估, 说明和用法请见微软帮助页面 PEARSON 函数. 如果要自己计算, 可以参考 wiki 的这个公式:

Pearson 相似度计算公式

C 的实现源码如下 (注意某些情况下可能会有计算精度丢失, 带来结果的不确定性)

double pearson_r(double x[], int x_n, double y[], int y_n) {
  if (x_n != y_n || x_n == 0) {
    return 0;
  }

  double n = (double)(x_n);
  double sum_x = 0;
  double sum_y = 0;
  double sum_x_sq = 0;
  double sum_y_sq = 0;
  double sum_x_by_y = 0;
  
  for (int i = 0; i < x_n; ++i) {
    sum_x += x[i];
    sum_y += y[i];
    sum_x_sq += x[i]*x[i];
    sum_y_sq += y[i]*y[i];
    sum_x_by_y += x[i]*y[i];
  }
  double res = n*sum_x_by_y - sum_x*sum_y;
  res /= sqrt(n*sum_x_sq - sum_x*sum_x);
  res /= sqrt(n*sum_y_sq - sum_y*sum_y);
  return res;
}