基于贝叶斯估计的星级得分排名
问题阐述
互联网早已成为人们生活的一部分，没事在网上看看电影、逛逛淘宝、定定外卖（有时间还是要多出去走走）。互联网的确为我们提供了非常多的便利，但它毕竟是一个虚拟的环境，具有更多的不确定性，大多数情况下我们只能通过别人的评论及打分来判别某个商品的好坏。五星打分是许多网站采用的商品排名方法，它也是消费者最直观最简单的评价尺度，我想大部分人都会去点击那些星级排名比较高的商品以最大限度降低我们的顾虑。
多数情况下，星级排名都能准确的反映一个商品的好坏，因为它是多人的一个综合得分，减小了个人偏好的影响。但是这里有一个前提条件，即打分的人要足够多。考虑下面两种商品：
1.A商品的星级平均得分为5（1位评论者）
2.B商品的星级平均得分为4.1（87位评论者）
A、B两种商品谁的得分排名更高呢？我想大部分都认为B应该排在A的前面吧，尽管B的平均星级得分要低于A，但是它有更多的体验人数，其得分更具有说服力。
现在我们已经明白，一个商品的排名应同时考虑它的星级得分与评论人数。那么我们应该如何将二者结合起来呢？也许会有人想到，我们可以给评论人数设定一个阀值，使得小于该阀值的商品，其排名会相对较低。上述过程可以使用下面的式子来表达：
这里的m代表平均星级得分，n代表打分的人数，k代表修正的阀值。K值该如何确定呢，它在某些极端的情况下准吗？这些都有待进一步考证。这里我们不采用该方法，我们运用贝叶斯大法。我们可以事先假设该商品事先满足某种分布（先验概率），并通过引入其实际打分情况，计算最终排名（后验概率）。
电影评论数据
本博文同样选取电影评分数据作为分析数据，它包含1664个电影，943个评分者，100000条打分记录。原始数据由明尼苏达大学社会计算研究中心提供，当然也可以从我的网盘下载已经处理好的数据。我们将下载后的数据放入data文件夹下。
为了管理的方便，本文分析的所有代码都将放入一个模块中，给该模块命名为stars.py。在该模块中，编写了一个名为Ratings的类，并把相关分析过程封装在该类中。
import os import numpy as np import pandas as pd import matplotlib.pyplot as plt from collections import Counter RATINGS = os.path.join(os.path.dirname(__file__), 'data', 'ratings.csv') class Ratings(): def __init__(self, path=RATINGS): self.path = path self.load() def load(self): self.data = pd.read_csv(self.path) def __str__(self): return str(self.data.head()) if __name__ == "__main__": ratings = Ratings() print ratings
上述打印结果表明，数据主要包含四个字段：user_id为客户的id号；movie_id为电影的id号，他与title是一一对应关系；rating为电影星级评分数据；title为电影的名称。
数据初探
数据导入后，我们就可以进行一些简单的分析了。依据电影的名称将所有电影进行分类，并计算每部电影的平均得分，取排名前十的电影进行展示。
class Ratings(object): . . . @property def movies(self): return self.data.groupby('title') def get_means(self): return self.movies['rating'].mean() def get_counts(self): return self.movies['rating'].count() def top_movies(self, n=-10): grid = pd.DataFrame({ 'mean': self.get_means(), 'conut': self.get_counts() }) return grid.ix[grid['mean'].argsort()[n:]] if __name__ == '__main__': ratings = Ratings() print ratings.top_movies()
上述结果是依据均值进行的电影排名，我们发现很多评分者较少的电影排在了非常靠前的位置，这显然不是我们想要的排名。
为了对均值评分的结果有一个更加直观的认识，我们使用密度分布图将其展现出来。
class Ratings(object): ... def plot_mean_frequency(self): grid = pd.DataFrame({ 'Mean Rating': self.movies['rating'].mean(), 'Number of Reviewers': self.movies['rating'].count() }) grid.plot(x='Number of Reviewers', y='Mean Rating', kind='hexbin', xscale='log', cmap='BuGn', gridsize=12, mincnt=1, title="Star Ratings by Simple Mean") plt.show()
上图结果表明，总体上所有电影的平均打分为要稍高于3.0。特别的，当打分人数非常少时，电影平均得分主要集中在4.0，3.0、1.0这几个值。当打分人数少于50时，平均得分为3.0左右。当打分人数大于50时，平均得分为3.5左右。
这里我们需要特别注意的是，当某个电影的评论人数非常少时，其平均得分具有较大波动性。这是由于当评论人数较少时，平均得分更容易受个别人的主观偏好的影响。而当评论人数大于100以后，电影平均得分主要集中在1.5至4.5之间，此时的评分结果更具有客观性。
从而，我们也不难看出，评分人数越多其平均得分就越高，这与我们实际情况也是相符的。在利用贝叶斯方法计算电影的后验得分时，我们可以假设其先验得分为整体的均值。
贝叶斯估计
贝叶斯估计主要是通过的引入最新的观察结果，不断修正我们的先验概率，从而求得最终的后验概率。我们首先来看一个简单贝叶斯的实例。假设我们要评估某个电影好评率（只考虑好坏两种情况），我们可以事先认为它的好评比例为0.5，即假设先验概率X为中性。然后我们在电影门口咨询N位观看者，统计他们各自的评论结果oi。最后我们就可以根据统计的结果O，不断修正好评率X的值，即求得X的后验概率。将上述过程用贝叶斯公式表示如下：
P(O)为观察结果O的概率分布值，它可以根据不同的X表现值，利用全概率公式求得。在不同X表现值下，P(O)是不变的，所以我们可以不做分析。P(X)为获取观察结果O之前，为X预先设定的概率分布，如果我们完全不了解X的分布情况，一般将其设定为均匀分布，我们也可以认为它为常值。现在我们将上述公式简化为如下形式：
回到前面的例子，假设好评概率为X，统计的N个电影评论结果中，有K个给予了好评，N-k个给予了差评，可以推导出如下后验概率的计算公式：
是不是发现上述公式变成了我们常见的伯努利分布。对于不同的N、K值，便可求得X的后验概率值，其分布图等价（这里只是意义上的等价，并不是它真正的概率密度分布）如下：
# -*- coding:utf-8 -*- from __future__ import division import numpy as np import matplotlib.pyplot as plt from functools import partial def f_NKX(N, K, X): return ((X)**K)*((1-X)**(N-K)) # 求X的偏函数 def f_NK(N, K): return partial(f_NKX, N, K) x = np.arange(0, 1, 0.01) # 分别设定不同的N、K值 p_list = [(1, 1), (5, 4), (10, 7), (25, 18), (50, 37), (100, 72)] plt.figure(facecolor='white') for i, p in enumerate(p_list): ax = plt.subplot('23'+str(i+1)) f = f_NK(p[0], p[1]) y = map(f, x) if i > 2: ax.set_xlabel('Probability of Prior') if i == 0 or i == 3: ax.set_ylabel('Probability of Posterior') ax.plot(x, y) plt.title('N=%d, k=%d' %(p[0], p[1])) plt.show()
从上图可以看出，随着观察数量N的增加，图形变得越来越集中，说明其后验概率期望值就会越来越明确，并随着观察数量结果的不断增加，后验期望概率就会无限接近实际观察的结果。这也再一次表明当评论人数足够多时，其结果的可信度将会更高，当评论人较少时其可信度将会较低。
Dirichlet分布
下面部分涉及的数学理论会比较多，倘若对Dirichlet分布不熟悉的建议先看《LDA数学八卦》。
在前面二项分布的例子中，并假设它的先验概率密度函数为均匀分布的情况。现考虑更一般的情况：假设其先验概率为Beta分布。为什么说Beta分布更具有一般性呢，先看下图：
百变星君Beta分布
观察上面的Beta分布图形，发现它是一个百变星君，可以是凹的、凸的、单调上升的、单调下降的；也可以是曲线，也可以是直线，而均匀分布也是特殊Beta分布的一种（图中的Beta(p|1,1)恰好就是就是均匀分布Uniform(0,1)）。正是由于Beta分布是如此的包罗万象，因此它在统计模型数据拟合和贝叶斯中应用广泛。
下面给出一般意义上的Beta分布函数（有没有发现它和Binomial分布很像）：
现在我们来了解Beat分布的一个非常重要的特性：Beta-Binomial共轭（具有相同形式的分布函数），《LDA数学八卦》是用如下式子表示的：
这里的加法形式容易给人造成误解，回到最原始贝叶斯估计的公式，发现先验概率与数据观察概率相乘才与后验分布是正相关的关系：
Beta-Binomial共轭有一个非常好的特性，当观察数据为二项分布时，参数的先验分布和后验分布都能保持Beta分布。
上述讨论的都是只是二项分布的情况，对于多项分布怎么办？前面提到的Dirichlet分布就派上用场了，一般形式的Dirichlet分布定义如下：
Dirichlet分布是Beta分布在多维上的推广。这里我也顺便给出多项分布Multinomial的公式：
和Beta-Binomial共轭一样，Dirichlet分布与Multinomial也存在共轭的情况：
Dirichlet分布还有一个非常方便的性质，其期望可以用下式求得：
回到电影星级评分
上一部分我们讨论了贝叶斯先验概率更一般的形式以及后验概率的估计过程，现在回到电影五星打分排名的问题。前面得知电影的得分可以是1,2,3,4,5星五种情况，现在我们的问题就是根据他们在五个星级水平下的得分比例，估计星级得分的后验概率分布，并求出最终期望得分。假设各个星级水平下的得分比例分别为P1、P2、P3、P4、P5，显然满足：
假设某部电影的N个评论者中打1,2,3,4,5星的人数分别为K1、K2、K3、K4、K5，由多项分布可得：
上式其实也是Dirichlet分布（X的表现形式是有限的情况下它与多项分布其实是等价的），其参数为：
这里我们假设P(X)的Dirichlet先验概率分布的系数为a0，那么后验概率可以表示为：
即满足如下参数的Dirichlet分布：
修正期望得分
求得后验概率分布的函数后，需要去解决最开始讨论的问题，即求得每一部电影最终的修正期望得分。由期望的相关知识可知E(wp)=wE(p)，w为常数，所以只要求得Dirichlet分布的期望概率就能很方便的求得期望得分，这里的权重即为星级得分水平。修正期望得分的求解表达式如下：
E(p)的求解公式在Dirichlet分布那部分讲解过了，这里就不在重复。
带入相关参数，修正期望得分排名得公式如下：
将上式化简可得：
仔细观察上式，会惊奇的发现，后验修正期望得分其实就是先验数据的总得分加上观察数据的总得分除以先验数据和观察数据之和。饶了一大圈原来我在推导调整平均公式（知道真相的我眼泪流下来），只不过这里的先验概率要事先“假设好”。现用更加贴近实际问题的方式来描述上述公式，设定参数m和c，其中：
- m代表先验分布的平均星级得分
- c代表先验分布中的评论人数
从而得到修正期望得分的贝叶斯估计为：
上述公式表明，当c越大时，该电影受先验分布的影响就会越大。从另一面也说明，该电影要想“摆脱”m的影响所需评论人数就会越多。上述过程的实现过程如下：
class Ratings(object): def __init__(self, path=RATINGS): self.path = path self.load() self.c = None self.m = None
. . . def set_prior_para(self, c, m): self.c = c self.m = m def bayesian_mean(self, arr): if not self.m or not self.c: raise TypeError("Bayesian mean must be computed with m and C") return (self.c * self.m + arr.sum()) / (self.c + arr.count()) def get_bayesian_estimates(self): return self.movies['rating'].agg(self.bayesian_mean) def top_movies(self, n=10): grid = pd.DataFrame({ 'mean': self.get_means(), 'count': self.get_counts(), 'bayes': self.get_bayesian_estimates(), }) return grid.ix[grid['bayes'].argsort()[-n:]] if __name__ == '__main__': ratings = Ratings() ratings.set_prior_para(10, 3) print ratings.top_movies()
我们用贝叶斯修正电影评分的主要目的是为了避免个人主观偏好的影响，而个人偏好在评论人数较多时会得到综合，所以c的值一般不要设的太大，设定先验参数c=10、m=3，观察其结果：
用贝叶斯修正的电影得分排名前十的结果表明，尽管某些电影的平均得分要高于其它电影的平均得分，但是由于它的观看人数较少，使得总排名靠后（观察第一行与第二行的结果）。
由于不同的先验分布的结果会导致不同排名结果，前面我们已经谈到c的值尽量不要太大，那么可不可以根据数据本身的特点来设定呢。我们可以筛选那些评论人数少于一定数量的电影，计算他们的平均评论人数与平均得分值，作为我们的先验参数c和m。在这里我将人数的阀值分别设定为5、10、50，观察其结果。
class Ratings(object):
. . .
def top_movies_by_num(self, n=10): grid = pd.DataFrame({ 'mean': self.get_means(), 'count': self.get_counts(), }) for i, num in enumerate([5, 10, 50]): filter_data = grid[grid['count'] <= num] # 计算评论人数少于num的所有电影的平均评论人数 self.c = np.mean(filter_data['count']) # 计算评论人数少于num的所有电影的平均得分 self.m = (filter_data['mean'].dot(filter_data['count']))/np.sum(filter_data['count']) grid['bayes%d'%(i+1)] = self.get_bayesian_estimates() return grid.ix[grid['bayes1'].argsort()[-n:]] if __name__ == '__main__': ratings = Ratings() print ratings.top_movies_by_num()
上述输出结果表明，选取不同阀值（阀值尽量不要过大）的观察数据作为贝叶斯先验分布的参数，其得分排序的结果会有所区别，但整体上都同时考虑了实际均值和评论人数的两个因素，一般的均值越大且评论人数越多的电影，其后验得分也会越大。
结论
本文讨论了电影评论得分的贝叶斯修正模型，并设定先验概率分布满足Dirichlet分布，有效的避免了个人偏好导致评分结果的片面性。使用电影五星打分数据进行实例分析，并使用python编程实现，其结果表明，基于贝叶斯修正的电影期望得分综合考虑了平均得分与评论人数，其排序结果更加符合实际结果。当然，本部分的工作还有待进一步完善，例如我们在打分时考虑某部电影的持续影响力等。虽然最后的结论比较简单，但还是花了我一些时间，希望能对大家有所帮助。
参考资料
