在广告系统中,一个重要的指标是CTR。ctr=点击(Click)/曝光(Impression)。
直觉上肯定是不对的,这种情况在微信大数据竞赛中也有出现,先来看一些数据。 【sum/count/ctr,图】
一个广告的每次展示,相当于概率试验里的投硬币。客户点击/不点击广告,是一个泊努利试验(Bernoulli Trial)。多个客户点击/不点击广告,可以看成是进行一串Bernoulli试验,即二项试验(Binomial Experiment),服从二项分布。
假设我们尝试了n次试验(有n个客户),每个客户点击该广告的概率为p(平均ctr)。该试验服从二项分布:B(n,p)。在n次试验中,观察到有k次点击的概率为:
例如,如果有100个visitors,该广告的点击率为10%,点击次数的概率分布(PMF)为:
即上面公式中:n=100, 横轴为k,纵轴为p。
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom
n, p = 4, 0.1
x = np.arange(1,100)
fig, ax = plt.subplots(1, 1)
ax.plot(x, binom.pmf(x, n, p), 'bo', ms=8, label='binom pmf')
plt.show()
我们可得到类似这样的曲线:
当CTR接近10%时,有相当高的机会看到8次-12次左右的点击。如果观察的次数越少,噪声也会越多。当只有4次观测(observation)时,有65%的机会看到没有点击,30%的机会看到一次点击,5%的机会看到2次点击。
from scipy.misc import comb
n = 4
k = 0
r = 0.1
p0 = comb(n,k)* (r)**k * (1-r)**(n-k)
### p1接近0.65
k = 1
p1 = comb(n,k)* (r)**k * (1-r)**(n-k)
### p2接近0.30
k=2
p2 = comb(n,k)* (r)**k * (1-r)**(n-k)
### p3接近0.05
假如,我们只有少量观测,是否可以去估计CTR呢?是否可以设计一个算法去模仿相应的模型数据?
为了在一个广告上模仿点击,我们首先使用一些分布上的CTR的值,接着使用它们作为在二项分布上的点击概率。这意味着我们需要两个随机变量。首先,是一个在[0,1]范围内的CTR的连续分布。第二个,是使用该CTR作为参数的二项分布。
绘制出对应的Beta曲线:
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
a,b=2,8
x = np.arange(0,1,0.01)
fig, ax = plt.subplots(1, 1)
ax.plot(x, beta.pdf(x, a, b), 'r-', lw=5, alpha=0.6, label='beta pdf')
plt.show()
可以得到类似的图:
从10次观测中的不确定性,通过该分布的扩展可以很容易解释。同理,我们可以尝试计算参数化,如果我们具有1000个观测以及200次点击。你可以注意到它将高度集中围绕在20%的CTR中。
将Beta分布与Binomial分布结合在一起,称为Beta-Binomial分布。
在参考文献一中,提出的方法是直接使用先验CTR:
通常,我们实际展示多个广告。计算观测时,当存在不确定性时,我们会生成一个CTR的估计值。一些广告可能具有成千上万的曝光和点击,而另一些可能只有很少的点击。这时候,可以引入贝叶斯推断(Bayesian inference)。我们使用一个先验的CTR,当新观测被记录时,我们不断更新它。先验CTR有很多方式确定。如果时间足够,我们可以使用基于Mean和sample size的参数化方法。如果我们没有足够信息,我们可以使用事先分配的先验CTR(non-informative prior),比如β(0.5,0.5):
一旦我们定义好先验CTR后,我们需要指定似然函数(likelihood function)。在我们的案例中,在结定参数集(CTR)下的观测的似然(likelihood)由二项分布给出。二项分布似然加上Beta先验,允许我们使用联合先验概率来获取一个后验分布。在我们的先验
下,经过N次曝光观测到x次点击,得到的后验是这样一个Beta分布:
。
还是使用之前的观测:4次曝光,1次点击,我们得到这样的后验:
贝叶斯模型可以后验的总结、或者通过随机抽样(均值、中位数、标准差等)进行分析。例如,在上面的贝叶斯推断步骤后,我们希望看到:一个广告在10000次曝光下,具有0.35的CPC(Cost-per-Click)
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta
cpc = 0.35
priorA = 0.5
priorB = 0.5
clicks = 1
impressions = 4
#rvs:生成随机数
plt.hist(beta.rvs(priorA + clicks, priorB + impressions - clicks, size=10000) * 10000 * cpc, 20, facecolor='green')
plt.xlabel('Cost')
plt.ylabel('Frequency')
plt.grid(True)
plt.show()
在许多情况下,我们更关心CTR的趋势,而非某个绝对值快照。对于罕见的(page/ad) pair,任何时间点,CTR的估计都是很大噪声的。根据某种数据连续性,我们可以将impression/click看做是重复测量的离散集合,并在这些测量上采用了一种指数平滑(exponential smoothing),最后我们得到一个平滑的CTR。
np.random.seed(0)
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 = np.random.beta(alpha, beta, num) #贝塔分布
I = []
C = []
for click_ratio in sample: #Beta分布生成的click_ratio
imp = random.random() * imp_upperbound
#imp = imp_upperbound
click = imp * click_ratio
I.append(imp)
C.append(click)
return I, C
def update_from_data_by_FPI(self, tries, success, iter_num, epsilon):
'''用不动点迭代法 更新Beta里的参数 alpha和beta'''
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 不动点迭代x_i+1 = g(x_i)'''
sumfenzialpha = 0.0
sumfenzibeta = 0.0
sumfenmu = 0.0
for i in range(len(tries)):
#special.digamma(z)是 在z处取gamma函数值 再求log
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)
def update_from_data_by_moment(self, tries, success): #tries是各组总样本数,success是各组click=1的样本和
'''用矩估计 更新Beta里的参数 alpha和beta'''
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):
'''矩估计'''
ctr_list = []
var = 0.0
for i in range(len(tries)):
ctr_list.append(float(success[i])/(tries[i] + 0.000000001))
mean = sum(ctr_list)/len(ctr_list)
for ctr in ctr_list:
var += pow(ctr-mean, 2) #方差
return mean, var/(len(ctr_list)-1)
具体调用代码:
HP = HyperParam(1, 1)
HP.update_from_data_by_moment(data['impression'].values, data['click'].values)
data['rate'] = (data['click'] + HP.alpha) / (data['impression'] + HP.alpha + HP.beta)