转化率预估中的贝叶斯平滑
1. 问题描述
在做比赛的过程中，我们发现了有转化率这个指标在大量数据下是有效的。理想情况下，例如某个广告点击量是10000次，转化量是100次，那转化率就是1%。但有时，例如某个广告点击量是2次，转化量是1次，这样算来转化率为50%。但此时这个指标在数学上是无效的。因为大数定律告诉我们，在试验不变的条件下，重复试验多次，随机事件的频率近似于它的概率。后者点击量只有2次，不满足“重复试验多次”的条件。
那么如何解决这个问题呢？
整体思路：用估计值来模拟实际CVR。
2. 解决方案
实际上，广告妆化率是随着时间迁移和用户喜好变化而变化的。因此我们可以利用先验知识来持续性地调整CVR。计算广告训练与平滑思想给出了一个很好的解决方案：贝叶斯平滑。
考虑到时序序列模型，我们把从第一天到第 天的所有先验知识汇总到第天的结果，并以此来对第 天的CTR进行平滑。在广告平滑上，没有什么方法比贝叶斯平滑能够更好的利用先验知识了，而帮助贝叶斯平滑方法实现目标的就是 分布。 分布的强大之处在于，通过改变其中的两个参数 和 ，你可以让 分布的图形变成任意形状，而且在加入先验知识前后，通过贝叶斯变换，对CTR的预估都可以表示为 分布。
分布中参数 和 的本质含义，即： 表示点击数， 表示曝光数。因为贝叶斯平滑的具体公式（后面再讲这个公式的原理）就是：
公式由来：
- 一般来说，点击还是不点击，这是服从伯努利二项分布的。
- 而二项分布的共轭分布就是 分布，也就是说，点击率服从 分布
- 我们可以从历史数据当中学到历史数据的 分布的具体参数，也就是先验分布 （不加任何条件限制的分布）
- 共轭先验有一个特性：如果找到一个 ，它是 的共轭先验，那么 的后验分布 和先验分布 会有一样的形式。
- 这个特性告诉我们：先验分布 （也就是历史数据）的分布与后验分布 （也就是 条件下点击率 的分布）是一个形式的
- 既然我们知道了先验分布 （也就是历史数据）的分布是 分布，那么我们就知道了后验分布 （也就是 条件下点击率 的分布）分布也是 分布
- 也就是说 ：先验分布 （也就是历史数据） + 后验知识 —> 后验分布 （也就是 条件下点击率的分布）
- 那么接下来我们就需要求解后验分布 的分布参数了
- 根据什么什么证明，后验分布的参数
- 其中I是展示数量，C是点击数量， 和 是历史数据的beta分布参数
- 那么后验分布\pi(r|x) 就是
- 如果我们要让点击率误差最小，那么取后验分布的均值，就是最好的点击率了！！！！也就是：
3. 参数估计
能不能直接把历史点击和历史曝光分别赋值给 和 来进行计算呢？显然不行，因为这么做就会犯之前我们提到的那些问题，比如不同日期的曝光、点击权重应该不一样。所以基础的贝叶斯平滑是不能解决我们刚才提到的问题的，我们需要深入研究 分布的特性，用一种新的方法通过先验知识求解 和 ，从而计算 。
参考文献: CTR预估中的贝叶斯平滑方法（二）参数估计和代码实现(https://www.bbsmax.com/A/A7zgmjRk54/)
3.1. 矩估计
矩估计的方法要追溯到19世纪的Karl Pearson，是基于一种简单的 “替换” 思想建立起来的一种估计方法。其基本思想是用样本矩估计总体矩. 由大数定理，如果未知参数和总体的某个(些)矩有关系，我们可以很自然地来构造未知参数的估计。具体计算步骤如下：
分布除了两个显性的重要参数 和 外，还有两个相对隐形但同样重要的参数，均值和方差，通过均值和方差可以唯一确定 和 的值，它们的数学关系如下：(见参考链接beta（贝塔）分布的数学期望和方差)
均值:
方差:
因此，如果我们根据数据集统计了平均值和方差，那么α和β的值也就确定了：
3.2. EM估计
这种估计方法其实叫Fixed-point iteration。只是有点类似EM的思想。
首先构造出似然函数，然后利用Fixed-point iteration来求得似然函数的最大值。
- 1）首先给出参数的一个初始值（通常可以使用矩估计得到的结果作为初始值）。
- 2）在初始值处，构造似然函数的一个紧的下界函数。这个下界函数可以求得其最大值处的闭式解，将此解作为新的估计用于下一次迭代中。
- 3）不断重复上述（2）的步骤，直至收敛。此时便可到达似然函数的stationary point。如果似然函数是convex的，那么此时就是唯一的最优解。
其实Fixed-point iteration的思想与EM类似。
#!/usr/bin/python
# coding=utf-8
import numpy
import random
import scipy.special as special
import pandas as pd
import time
import math
from math import log
class HyperParam(object):
def __init__(self, alpha, beta):
self.alpha = alpha
self.beta = beta
def sample_from_beta(self, alpha, beta, num, imp_upperbound):
sample = numpy.random.beta(alpha, beta, num)
I = []
C = []
for click_ratio in sample:
imp = random.random() * imp_upperbound
#imp = imp_upperbound
click = imp * click_ratio
I.append(imp)
C.append(click)
return I, C
# 平滑方式1
def update_from_data_by_FPI(self, tries, success, iter_num, epsilon):
'''estimate alpha, beta using fixed point iteration'''
'''tries ： 展示次数
success : 点击次数
iter_num : 迭代次数
epsilon : 精度
'''
for i in range(iter_num):
new_alpha, new_beta = self.__fixed_point_iteration(tries, success, self.alpha, self.beta)
# 当迭代稳定时，停止迭代
if abs(new_alpha-self.alpha)<epsilon and abs(new_beta-self.beta)<epsilon:
break
self.alpha = new_alpha
self.beta = new_beta
def __fixed_point_iteration(self, tries, success, alpha, beta):
'''fixed point iteration'''
sumfenzialpha = 0.0
sumfenzibeta = 0.0
sumfenmu = 0.0
# digamma 指伽马函数，是阶乘在实数与复数域的扩展
sumfenzialpha = special.digamma(success+alpha) - special.digamma(alpha)
print sumfenzialpha
# for i in range(len(tries)):
# sumfenzialpha += (special.digamma(success[i]+alpha) - special.digamma(alpha))
# sumfenzibeta += (special.digamma(tries[i]-success[i]+beta) - special.digamma(beta))
# sumfenmu += (special.digamma(tries[i]+alpha+beta) - special.digamma(alpha+beta))
return alpha*(sumfenzialpha/sumfenmu), beta*(sumfenzibeta/sumfenmu)
# 平滑方式2
def update_from_data_by_moment(self, tries, success):
'''estimate alpha, beta using moment estimation'''
# 求均值和方差
mean, var = self.__compute_moment(tries, success)
#print 'mean and variance: ', mean, var
#self.alpha = mean*(mean*(1-mean)/(var+0.000001)-1)
self.alpha = (mean+0.000001) * ((mean+0.000001) * (1.000001 - mean) / (var+0.000001) - 1)
#self.beta = (1-mean)*(mean*(1-mean)/(var+0.000001)-1)
self.beta = (1.000001 - mean) * ((mean+0.000001) * (1.000001 - mean) / (var+0.000001) - 1)
def __compute_moment(self, tries, success):
# 求均值和方差
'''moment estimation'''
ctr_list = []
# var = 0.0
mean = (success / tries).mean()
if len(tries) == 1:
var = 0
else:
var = (success / tries).var()
# for i in range(len(tries)):
# ctr_list.append(float(success[i])/tries[i])
# mean = sum(ctr_list)/len(ctr_list)
# for ctr in ctr_list:
# var += pow(ctr-mean, 2)
return mean, var
def test():
#设定初始值
hyper = HyperParam(1, 1)
#--------sample training data--------
# I, C = hyper.sample_from_beta(10, 1000, 10000, 1000)
# print I, C
train_data = pd.read_csv('data.csv',nrows=10000)
print 'read finish'
# 统计点击次数和转化次数
key = ['creativeID']
train_data['count'] = 1
train_data = train_data.groupby(key).agg('sum').reset_index()
# 此时，train_data['count']是点击次数
# train_data['label']是点击次数
print 'cal finish'
I = train_data['count']
C = train_data['label']
print key
start = time.clock()
#--------estimate parameter using fixed-point iteration--------
# 计算平滑
hyper.update_from_data_by_FPI(I, C, 1000, 0.00000001)
end = time.clock()
print hyper.alpha, hyper.beta
print 'run time: ',end - start
start1 = time.clock()
#--------estimate parameter using moment estimation--------
hyper.update_from_data_by_moment(I, C)
end1 = time.clock()
print hyper.alpha, hyper.beta
print 'EM run time: ', end1 - start1
if __name__ == '__main__':
test()4.贝叶斯估计
参考文献转化率（CTR）预测的贝叶斯平滑: https://blog.csdn.net/jinping_shi/article/details/78334362
4.1. 点击率平滑的假设
对于一件商品或一条广告，对于某次曝光，用户要么点击，要么没点击，这符合二项分布。因此下文中对于点击率类的贝叶斯平滑，都是基于以下假设：对于某件商品或广告XX，其是否被点击是一个伯努利分布（Bernoulli）。
其中 表示某个广告是否被点击，点击取1，未被点击取0， 是某件商品被点击的概率，即点击率。
4.2. 冷启动问题——点击率极大似然估计
在上式的假设下，可以使用极大似然法计算出点击率的估计 。具体做法为：
- 从用户日志中随机抽取 条记录，对任一条记录 都有
那么所有记录的点击数的联合概率密度就是上式的连乘，也就是构造了极大似然函数。将极大似然函数对 求导并令导数等于0，就可以解出 的估计值 。 就是点击率的极大似然估计。当某个商品的点击次数或者曝光次数为0时，可以用 当成它的初始值。
然而这样并没有解决新上线广告问题。例如有两件商品A和B，其点击率分别为和，，但商品的曝光只有10次，商品的曝光有100次，这样比较是否合理？那么我们就要用到贝叶斯估计来解决这个问题了！
4.3. 广告投放不足问题——点击率的贝叶斯估计
在贝叶斯框架下，我们假设点击率 服从某个分布：
因为这是基于经验的，这个分布称为先验分布。贝叶斯参数估计可以同时解决最开始提出的两个问题。其过程是基于经验或历史数据先给出一个的估计值，然后基于现有的数据在这个估计值上修正，得到最终的点击率估计，此时的估计值已经是修正过的。更美好的是我们可以分离出修正参数，来进行更好的估计（即3.1 公式中的和）。
既然有先验分布，就有后验分布。的后验分布记作 。其中表示输入数据或特征，对于点击率预测，就是点击次数和曝光量。因为已经看到了数据，才确定的分布，因此叫做『后验』分布。贝叶斯估计的实质就是求后验分布。即基于当前的点击次数和曝光量，求点击率的分布；而未看到数据之前点击率的分布是。下面会讲解如何计算后验分布.
贝叶斯估计的过程可以简单认为：
用历史数据根据估计，记作 ；用当前数据根据估计 ，记作 current，然后用 修正 。
4.4. 损失函数
的后验分布 是个概率密度函数，无法知道 确切的值。需要求出最接近真实情况的 需要损失函数来约束。
适用于点击率的损失函数有：
贝叶斯参数估计的过程可以简单描述为：求 ，使得损失函数在r的后验分布上的期望最小。这句话的数学定义为：
因此需要知道的形式，然而的形式一般不知道的，但是可以知道的形式（经验值嘛，我们指定的）。此外，数据的分布我们也是知道的，其概率密度函数（pdf）记为，表示数据的分布跟参数有关，是要求解的参数，在这里就是点击率。
这时可以根据贝叶斯公式计算出：
其中表示边缘概率密度定义，
上式好复杂，但其实一些常见的分布都可以求出上式积分的具体形式。但通常不用实际去积分，因为满足一定条件，跟有一样的形式。是已知的，如果形式一样，也就容易求得了。下面介绍共轭先验的概念。
共轭先验：如果找到一个，它是的共轭先验，那么的后验分布和先验分布会有一样的形式。
『轭』是指驾车时套在牲口脖子上的曲木。古代拉扯的牲口通常有两只，因此轭是连接两只牲口的工具。在这里共轭是指和通过联系起来了。
之前假设广告是否点击服从伯努利分布，参数为；对于点击次数，服从的是二项分布，即 ，二项分布的共轭先验是分布。分布的参数是和，即αβ ,根据共轭先验的定义，的后验分布的形式跟其先验分布一样，即αβ。
对于点击率预测，求出，带入公式4.4 中的 公式，当时，
上式的求解过程可以参考贝叶斯参数估计最后的例子。上式就是点击率估计（平滑）的最终形式。其中和就是点击次数和曝光量，即为3.2中的，αβ是3.2中的。和是从历史数据中得到的。
上面的内容给出了为什么很多文章会假设点击率服从分布的理由，因为最终的平滑的因子是分布（先验分布）中的两个参数。