机器学习手记

补遗: 机器学习手记系列 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

为什么要预估点击率

背景

想到这个题目是因为 @lijiefei 某天跟我说他有师弟面淘宝时被问到 “点击率预估的目标到底是什么”, 笨狗当时胡乱扯了一通, 发现要把这个似乎已经是真理的事情掰清楚还没那么容易, 于是有此念想写文一篇详细分析下原因

我和 jiefei 认识是在百度做搜索广告的时候, 那就从搜索广告开始说为什么要预估点击率, 以及预估点击率的目标. 先申明一些名词和假定:
1) 每个广告 (Ad) 有一个出价 (Bid), 并有其在某情形下实际的点击率 (Click-Through-Rate, CTR)
2) 广告按点击收费 (Charge per Click, CPC), 下面我们会分别讨论一价计费 (First-Price, FP, 即广告出价多少则一次点击计费多少) 和二价计费 (Second-Price, SP, 即广告按下一位出价来支付点击价格, 更普遍的是 GSP)
3) 千次展现收费 (Cost Per Mille, CPM, 或 RPM, R for Revenue), 即对点击付费广告其展示一千次情况下的收入 (一价计费下等价于 1000*CTR*Bid), 或是展示广告的千次展现固定价格
4) 预估点击率 (predict CTR, pCTR) 是指对某个广告将要在某个情形下展现前, 系统预估其可能的点击概率

目标分类

搜索广告跟自然结果一个很大的区别就是自然结果只要有一点相关就应该放到所有结果里去, 至于先后位置那个再说, 而广告, 是有个相关性的准入门槛的, 不相关的广告出价再高, 丢出来还是会被骂死. 那怎么判断相关? 用户会用鼠标点击来对结果投票, 相关的广告会被点击, 不相关的广告不会被点击, 那很自然就能得出 “点击率和相关性正相关” 这个结论 (至于描述里写 “二十五岁以下免进” 但实际是钢材广告的这种诱骗行为后面再说怎么处理). 那对于这种相关性准入的场景, 预估点击率就是预估广告是否相关, 最朴素情况下这是个二分类问题, 那不管预估成怎样, 只要有一种分割方法能分开是否相关就行了. 此时预估点击率的目标是能对广告按相关与否分类 (或说按相关性排序并给出一个截断值). 评估分类问题好坏, 一般都是看准确和召回两个指标, 用人工打分的记录来做回归验证就行

目标排序

判断相关与否只是点击率预估对广告的一个小辅助, 我们来看看广告的目标是什么? 没错, 是赚钱. (我曾经在其他场合说过广告的目标是维持用户体验下持续赚钱, 不过跟赚钱这一简化目标这不冲突, 前面相关性上已经保证了维持用户体验, 那只要能让广告主还有的赚, 就能持续赚钱) 我们再把问题简化下, 如果广告都是一样的固定价格, 且就以这个价格按点击计费, 那在 PV 一定且预算充分的情况下, 更高的点击率则意味着更赚钱. 这样目标可以等价于怎么挑出更赚钱的广告, 就是那些点击率最高的广告, 我们只要能弄明白广告实际点击率的高低关系就能取得收益最大化, 预估点击率在这时候又是个排序问题, 我们只要弄对广告之间的序关系, 就可以收益最大. 评估排序问题的好坏, 一个经典方法是对 pCTR 的 ROC 曲线算 AUC (曲线下面积), 实际上我见过的做法也都是通过评估 AUC 的高低来判断点击率预估模型的好坏

目标带权排序

上一段里对广告这个业务做了很多简化, 比如大家价格都是一样的, 如果我们考虑价格不一样的情况, 那预期收益就会变成 (价格Bid*点击率CTR), 这个值很多地方也叫 CPM 或 RPM. 如果是对 CPM 排序, 那就需要我们预估的点击率在维持序关系正确的前提下, 还要保证相互之间的缩放比是一样的. 比如有广告 A, B, C, 实际点击率是 5%, 3%, 1%, 那在价格一致的情况下, 我预估成 5-3-1 还是 5-4-3 是没关系的, 但在价格不一样的情况下, 比如 1, 1.5, 3, 这时候 5-4-3 的预估点击率值会让他们的预估排名和实际排名刚好颠倒过来, 不过预估 5-3-1 或 10-6-2 (放大一倍) 倒没关系. 为了评估这个结果, 可以在描 ROC 曲线时把价格乘上去, 那最后还是判断排序问题的好坏, 加了价格的 AUC 我们可以叫 wAUC (weighted-AUC), 这个离线评估和在线效果依然可以对等

目标准确

从准确召回到 AUC 再到 wAUC, 看起来对已有问题可以完美解决了? 不过广告计费显然不是 FP 这么简单, 在 Google 的带领下几乎所有的搜索引擎都使用了 GSP (Generalized Second Price) 来对广告点击进行计费, 这里再简单解释下最简版 GSP 是怎么回事:
1) 所有广告按 CPM 排序 (即 CTR*Bid)
2) 排最后的广告收底价 (Reserved Price, RP)
3) 其他广告按他的下一位 CPM_i+1 除以他的 CTR_i 并加一个偏移量来计费, 并保证比底价高, 即 Price_i = max(CPM_i+1/CTR_i + delta_price, RP)

至于 GSP 的细节和为什么这么做能保证收入和体验的平衡等可以详见相关论文, 我们只讨论在 GSP 模式下, 点击率预估的作用和关键点. 根据介绍 GSP 时最后那个公式, 如果把 CPM_i+1 拆成 CTR_i+1*Bid_i+1, 看起来只要保证同比缩放还是会没问题? 但是, 凡事怕但是, 在搜索广告里, 不同的展现位置对点击率还有影响, 比如广告 A, B 在第一位点击率是 5%, 3%, 而在第二位是 3%, 2%, 那只是同比缩放就很难保证最终比较是一致的问题了, 所以最好还是保证预估值跟实际值尽可能接近的好, 这样才能在预估时获得更实际用时完全一样的场景. 评估准确度, 我们有 MAE 和 MSE 等一堆指标, 也是现成的工作的比较好的东西

扩展和吐槽

有行家可能会吐槽说我刚那个不同广告在不同位置的衰减不一致这个说法, 跟公开论文说的不一样, Yahoo 的 paper 里说不同广告在同位置的衰减是一样的. 我只能说, 骚年, 你太天真了… 衰减因子怎么可能只是 f(pos) 这样一个简单函数, 从实际情况来看, 衰减函数和广告是有关的, 但我们又不能对每个广告都去估一个 f(pos, ad), 好在, 我们发现可以把不同的广告做聚类后得到一个 f(pos, type) 的函数簇, 事实上, 最后的衰减函数不仅仅有 pos 和 type 两个因子, 而且里面的因子可以极度简化, 最后的衰减用简单函数就能很好拟合, 我说的够多了, 再说估计要被前东家找麻烦, 你们来感受一下就好

前面也提到介绍的 GSP 是最简版的, 那如果是正常版会有什么不一样? 那就是排序时用的是 Quality*Bid, 这个 Quality 百度叫质量度, 也就是广大广告主望眼欲穿的那几颗星. Quality 和 CTR 有什么不一样? 最简情况下这两个当然是一样的, 但是我们可能还要考虑广告主的信誉度, 目标网页的质量等因素 (比如前面提到过的那种描述欺诈或诱导的广告在这个 Quality 因子里被调整掉), 最后这个 Quality 就会是包含了 CTR 的一个多因子复合值, 那要准确估计这个复合值, 当然也要求其中的每个因子的值尽可能准确. 这里在炒冷饭说准确度, 以及 MAE/MSE 的作用

实际上据我所知各搜索广告平台用的是比正常版的 GSP 还加强或改造过的版本, 里面的因子, 公式, 逻辑更复杂. 在这种情况下还是需要继续强调 CTR 预估的准, 才能做更精确的预估, 从而带来更大的收益

广告之外

呼~ 说完了搜索广告, 我们再简单说下内容广告. 搜索广告几乎都是点击付费, 而内容广告同时存在按点击付费和按展现付费, 那怎么比较一个按点击付费的广告和一个按展现付费的广告哪个预期收益更高? 同样是 CPM, 按展现付费的广告给的是确定值, 而按点击付费的是一个预估值, 通过 Bid*pCTR 得到, 如果 pCTR 不准, 就会导致点击付费广告的预期收益计算不准, 则不一定受益最大. 继续强调预估的准的好处

说完广告我们就能说说其他的, 比如搜索, 比如推荐, 这几个的优化目标如果是带来的量, 因为总体 PV 我们没法人工干预, 且每个点击是等价的, 那最后的优化就变成了点击率, 预估的序关系越接近真实, 可获得的收益更高. 如果不同的点击价值不一样, 那就可以把这个点击换做价格代入广告的模型, 因为没有二价计费那么讨厌的变换, 所以就按一价计费来考虑, 保证序正确且等比缩放就能保证收益最大. 如果再激进一点, 评估收益时还加入更复杂的因子而不仅仅是价格这个独立因素, 那自然就要求点击率预估准确, 从而保证做决策时和实际情况一致, 继而保证最终的优化目标最大化

总结

1) 点击率预估是为产品的最终目标服务的, 最终目标可以是广告的收入, 广告的相关性, 推荐的接受率等, 看具体场景
2) 点击率预估的直接目标根据需求场景不同, 分别是保证预估值和实际值分类正确, 预估序和实际序正确, 预估值和实际值是等比缩放的, 预估值等于实际值
3) 要保证离线评估点击率预估的效果, 分别可用分类的准确率和召回率, 排序的 AUC, 带权排序的 wAUC, 相似度 MAE/MSE 来评估

机器学习手记系列 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;
}