本文用到的包
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats
import math
TrueSkill算法是Elo排名方法与贝叶斯规则的结合，可用于计算竞赛选手的能力排名。文献1提出了这个方法，文献2创造性地建议使用这方法来计算问答类社区问题的难度。
算法给每一个用户分配一个正态分布，均值代表真实能力，方差代表系统对该用户真实能力猜测的不确定程度。一开始假设每个人分布的均值和方差一致，此后利用数据中的每一对输赢关系不断更新每个用户的分布。更新规则如下：
图1. TrueSkill规则一
其中
图2. TrueSkill规则二
方程（7）中的分子N函数和分母Phi函数，分别是标准正态分布的PDF与CDF。
算法的精髓是消灭数据中的惊奇（或者充分吸收数据信息）。如果系统中储存的分布是A的均值大于B的，数据中却出现了B赢了A的情况，那么就把B的均值大幅提升，把A的均值大幅减小；如果不出意外地，A赢了B，那么把A的均值小幅提升，把B的均值小幅减小。
图3. TrueSkill规则展示
上图中蓝色和棕色分布分别对应A和B两个用户，实线是原来的分布，虚线是根据一个新的输赢关系更新后的分布。左下角图展示了A赢了B的情况，右下图展示B赢了A的情况。上面两个图展示了，输赢关系带来均值（左图）和方差（右图）的改变的剧烈程度。t<0的时候是令人惊奇的，因为此时系统中储存的赢家的均值比输家还要低。越惊奇，改变越大。
该算法的代码：
# updataing function with a beta parameter, bigger beta updates more dramatically
v = lambda x: scipy.stats.norm(0,1).pdf(x)/scipy.stats.norm(0,1).cdf(x)
w = lambda y: v(y)*(v(y) + y)
def update(mw,sw,ml,sl,beta): # mw: miu of winner; sw: sigma of winner; ml: miu of loser; sl: sigma of loser;
t = mw - ml
c = np.sqrt(2*beta**2+sw**2+sl**2)
mw += v(t/c)*sw**2/c
ml -= v(t/c)*sl**2/c
sw *= np.sqrt(1 - w(t/c)*sw**2/c**2)
sl *= np.sqrt(1 - w(t/c)*sl**2/c**2)
return mw,sw,ml,sl
绘制图3的代码如下
准备工作
# normal distribution ploting function
def plotNormalPDFs(xmin,xmax,m1,s1,m2,s2,col1,col2,line,ax,lab1,lab2):
x = np.linspace(xmin,xmax,200)
y1 = map(lambda x:scipy.stats.norm(m1, s1).pdf(x), x)
y2 = map(lambda x:scipy.stats.norm(m2, s2).pdf(x), x)
ax.plot(x,y1,linestyle=line,color=col1,label=lab1,linewidth=1.5)
ax.plot(x,y2,linestyle=line,color=col2,label=lab2,linewidth=1.5)
# generate data
xs=np.linspace(-6,6,100)
vs=map(lambda x:v(x),xs)
ws=map(lambda x:w(x),xs)
xmin,xmax,m1,s1,m2,s2 = 1,100,60,5,40,20
beta = (s1+s2)/4
m1w,s1w,m2l,s2l = update(m1,s1,m2,s2,beta)
m2w,s2w,m1l,s1l = update(m2,s2,m1,s1,beta)
画图代码如下。因为简书编辑器对代码中的$号敏感，会自动识别为复制段落，所以这里去除了横轴和纵轴标签中latex语法的$号，实际绘制中需补上，例如将
ax1.set_ylabel(r'v(t)',size=14)
修订为ax1.set_ylabel(r'$v(t)$',size=14)
# plot
fig = plt.figure(figsize=(10, 8),facecolor='white')
ax1 = fig.add_subplot(2,2,1)
ax1.plot(xs,vs, color='YellowGreen',linewidth=1.5)
ax1.set_xlabel(r't = \mu_{winner} - \mu_{loser}',size=14)
ax1.set_ylabel(r'v(t)',size=14)
#
ax2 = fig.add_subplot(2,2,2)
ax2.plot(xs,ws, color='Tomato',linewidth=1.5)
ax2.set_xlabel(r't = \mu_{winner} - \mu_{loser}',size=14)
ax2.set_ylabel(r'w(t)',size=14)
#
ax3 = fig.add_subplot(2,2,3)
col1 = 'SteelBlue'; col2='Chocolate'; lineo='-'; linen='-.'
plotNormalPDFs(xmin,xmax,m1,s1,m2,s2,col1,col2,lineo,ax3,'original i ','original j')
plotNormalPDFs(xmin,xmax,m1w,s1w,m2l,s2l,col1,col2,linen,ax3,'updated i','updated j')
ax3.set_xlabel(r'skills',size=14)
ax3.set_ylabel(r'propability',size=14)
ax3.legend(fontsize=10,loc=2)
#
ax4 = fig.add_subplot(2,2,4)
plotNormalPDFs(xmin,xmax,m1,s1,m2,s2,col1,col2,lineo,ax4,'original i ','original j')
plotNormalPDFs(xmin,xmax,m2w,s2w,m1l,s1l,col2,col1,linen,ax4,'updated j','updated i')
ax4.set_xlabel(r'skills',size=14)
ax4.set_ylabel(r'propability',size=14)
ax4.legend(fontsize=10,loc=2)
#
plt.tight_layout()
plt.show()
#plt.savefig('xxx/demo.png',transparent=True)
最后，如果进行大规模计算，可以使用如下版本的update函数。其对正态分布函数预先进行了数值计算，不用每次都生成新的正态分布函数，因此大大提升速度
# fast version of updating function
def fastupdate(mw,sw,ml,sl,beta): # miu and sigma of winner and loser
sw2=math.pow(sw,2)
sl2=math.pow(sl,2)
t = mw - ml
c = math.sqrt(34.72222+sw2+sl2)
c2=math.pow(c,2)
tc = t/c
vtc = 0.79788*0.60653**math.pow(tc,2)/math.erfc(-0.70711*tc)
wtc = vtc*(vtc + tc)
mw += vtc*sw2/c
ml -= vtc*sl2/c
sw *= math.sqrt(1 - wtc*sw2/c2)
sl *= math.sqrt(1 - wtc*sl2/c2)
return mw,sw,ml,sl
