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。若是轉載，請註明出處，在未經做者贊成下將本文用於商業用途，將追究其法律責任。