CTR预估中的贝叶斯平滑方法（二）参数估计和代码实现

1. 前言

前面博客介绍了CTR预估中的贝叶斯平滑方法的原理http://www.cnblogs.com/bentuwuying/p/6389222.html。

这篇博客主要是介绍如何对贝叶斯平滑的参数进行估计，以及具体的代码实现。

首先，咱们回顾一下前文中介绍的似然函数，也就是咱们须要进行最大化的目标函数：

下面咱们就基于这个目标函数介绍怎样估计参数。

 

2. 参数估计的几种方法

1. 矩估计

矩估计在这里有点乱入的意思：），由于它其实不是用来最大化似然函数的，而是直接进行参数的近似估计。

矩估计的方法要追溯到19世纪的Karl Pearson，是基于一种简单的 “替换” 思想创建起来的一种估计方法。 其基本思想是用样本矩估计整体矩. 由大数定理，若是未知参数和整体的某个(些)矩有关系，咱们能够很天然地来构造未知参数的估计。具体计算步骤以下：

1）根据给出的几率密度函数，计算整体的原点矩（若是只有一个参数只要计算一阶原点矩，若是有两个参数要计算一阶和二阶）。因为有参数这里获得的都是带有参数的式子。好比，有两个参数时，须要先计算出：指望  ; 方差  。在Beta分布中，能够计算获得，E(x) = α / (α+β)，D(x) = αβ / (α+β)²(α+β+1)。

2）根据给出的样本，按照计算样本的原点矩。一般它的均值mean用  表示，方差var用  表示。（另外提一句，求时，一般用n-1为底。这样是想让结果跟接近整体的方差，又称为无偏估计）

3）让整体的原点矩与样本的原点矩相等，解出参数。所得结果即为参数的矩估计值。这里有，mean = E(x) = α / (α+β)，var = D(x) = αβ / (α+β)²(α+β+1)。因而乎，咱们能够求得α，β：

α = [mean*(1-mean)/var - 1] * mean

β = [mean*(1-mean)/var - 1] * (1-mean)

 

2. Fixed-point iteration

首先构造出似然函数，而后利用Fixed-point iteration来求得似然函数的最大值。

1）首先给出参数的一个初始值（一般能够使用矩估计获得的结果做为初始值）。

2）在初始值处，构造似然函数的一个紧的下界函数。这个下界函数能够求得其最大值处的闭式解，将此解做为新的估计用于下一次迭代中。

3）不断重复上述（2）的步骤，直至收敛。此时即可到达似然函数的stationary point。若是似然函数是convex的，那么此时就是惟一的最优解。

其实Fixed-point iteration的思想与EM相似。

首先给出两个不等式关系：

1）

2）

由此能够获得对数似然函数的一个下界：

 

想要获得此下界函数的最大值，能够分别对α，β求偏导，并令之等于零，此时便获得α和β各自的迭代公式：

 由此，每次迭代，参数都会达到这次下界函数的最大值处，同时也就使得对应的似然函数值也相应地不断增大，直至收敛到似然函数的最大值处。

 

3. EM

经过将几率参数做为隐含变量，任何估计几率参数的算法均可以使用EM进一步变成估计个数参数的算法。

（1）E-step：计算出隐含变量p在已观测数据（观测到的每一个类别发生的次数，以及每一个类别的超参数值的上一轮迭代的取值）下的后验分布，即可以获得彻底数据的对数似然函数的指望值。

（2）M-step：对E-step中的指望值求最大值，即可获得相应的超参数的本轮迭代的更新值。

（3）不断重复地运行E-step和M-step，直至收敛。

回来咱们这里的问题上，在有了观测到的每一个类别发生的次数，以及每一个类别的超参数值的上一轮迭代的取值后，隐含变量p的后验分布为：

而此时的彻底数据的对数似然函数的指望为：

其中，

因而乎，咱们能够对彻底数据的对数似然函数的指望求最大值，从而获得α，β的更新值，有不少方法，直接求偏导，梯度降低，牛顿法等。

可是呢，此时咱们并不须要很是精确地求得它的最大值，而是仅仅用牛顿法迭代一次。相比于精确地求得最大值，这种方法在每次迭代时只有一半的计算量，可是迭代次数会超过两倍。

牛顿法的迭代可见：

 

3. 代码实现

 1 #!/usr/bin/python
 2 # coding=utf-8
 3 
 4 import numpy
 5 import random
 6 import scipy.special as special
 7 import math
 8 from math import log
 9 
10 
11 class HyperParam(object):
12     def __init__(self, alpha, beta):
13         self.alpha = alpha
14         self.beta = beta
15 
16     def sample_from_beta(self, alpha, beta, num, imp_upperbound):
17         sample = numpy.random.beta(alpha, beta, num)
18         I = []
19         C = []
20         for click_ratio in sample:
21             imp = random.random() * imp_upperbound
22             #imp = imp_upperbound
23             click = imp * click_ratio
24             I.append(imp)
25             C.append(click)
26         return I, C
27 
28     def update_from_data_by_FPI(self, tries, success, iter_num, epsilon):
29         '''estimate alpha, beta using fixed point iteration'''
30         for i in range(iter_num):
31             new_alpha, new_beta = self.__fixed_point_iteration(tries, success, self.alpha, self.beta)
32             if abs(new_alpha-self.alpha)<epsilon and abs(new_beta-self.beta)<epsilon:
33                 break
34             self.alpha = new_alpha
35             self.beta = new_beta
36 
37     def __fixed_point_iteration(self, tries, success, alpha, beta):
38         '''fixed point iteration'''
39         sumfenzialpha = 0.0
40         sumfenzibeta = 0.0
41         sumfenmu = 0.0
42         for i in range(len(tries)):
43             sumfenzialpha += (special.digamma(success[i]+alpha) - special.digamma(alpha))
44             sumfenzibeta += (special.digamma(tries[i]-success[i]+beta) - special.digamma(beta))
45             sumfenmu += (special.digamma(tries[i]+alpha+beta) - special.digamma(alpha+beta))
46 
47         return alpha*(sumfenzialpha/sumfenmu), beta*(sumfenzibeta/sumfenmu)
48 
49     def update_from_data_by_moment(self, tries, success):
50         '''estimate alpha, beta using moment estimation'''
51         mean, var = self.__compute_moment(tries, success)
52         #print 'mean and variance: ', mean, var
53         #self.alpha = mean*(mean*(1-mean)/(var+0.000001)-1)
54         self.alpha = (mean+0.000001) * ((mean+0.000001) * (1.000001 - mean) / (var+0.000001) - 1)
55         #self.beta = (1-mean)*(mean*(1-mean)/(var+0.000001)-1)
56         self.beta = (1.000001 - mean) * ((mean+0.000001) * (1.000001 - mean) / (var+0.000001) - 1)
57 
58     def __compute_moment(self, tries, success):
59         '''moment estimation'''
60         ctr_list = []
61         var = 0.0
62         for i in range(len(tries)):
63             ctr_list.append(float(success[i])/tries[i])
64         mean = sum(ctr_list)/len(ctr_list)
65         for ctr in ctr_list:
66             var += pow(ctr-mean, 2)
67 
68         return mean, var/(len(ctr_list)-1)
69 
70 
71 
72 def test():
73     hyper = HyperParam(1, 1)
74     #--------sample training data--------
75     I, C = hyper.sample_from_beta(10, 1000, 10000, 1000)
76     print I, C
77 
78     #--------estimate parameter using fixed-point iteration--------
79     hyper.update_from_data_by_FPI(I, C, 1000, 0.00000001)
80     print hyper.alpha, hyper.beta
81 
82     #--------estimate parameter using moment estimation--------
83     hyper.update_from_data_by_moment(I, C)
84     print hyper.alpha, hyper.beta

 

4. 参考文献

1. Click-Through Rate Estimation for Rare Events in Online Advertising

2. Estimating a Dirichlet distribution

 

版权声明：

   本文由笨兔勿应全部，发布于http://www.cnblogs.com/bentuwuying。若是转载，请注明出处，在未经做者赞成下将本文用于商业用途，将追究其法律责任。