机器学习项目中不可忽视的一个密辛 - 大数定理

1. 大数定理是什么?

大数定理是一个非常底层的基础性原理,大量的机器学习理论和算法实际上都建立在其之上。只是因为其太底层,也太合理的,合理到我们感觉这就是一个毫无争议的真理而不需要去单独讨论。我们本文和读者来一起从贝叶斯推断的角度来尝试对大数定理展开一些讨论。

0x1:大数定理基本数学定义

设𝑍𝑖 为服从特定概率分布的𝑁个独立采样样本,按照大数定理理论,只要𝑍的期望𝐸[𝑍]是有限的,则就有如下关系:

用文字来描述:即来自同一分布的一组随机变量,其均值收敛域该分布的期望

以上等式只有在极限条件下成立,但是随着参加平均计算的样本越来越多,均值将越来越接近𝑍的期望值𝐸[𝑍]。除了一些极少数的特殊情况外,大数定理对于大多数分布都是正确的。

0x2:泊松随机变量的收敛 - 大数定理普遍存在的一个例证

下面是三个不同的泊松随机变量在大数定理上的实例,用三个不同的曲线表示。 三个泊松分布的参数𝜆 = 4.5,样本个数 sample_size=100000 个。

泊松分布本身也是一个随机过程,参数 𝜆 一样,只是表明它们符合同一个概率分布,但是泊松分布本身每次的取值依然是一个随机过程。这里用三个泊松分布来进行实验,只是想说明不管过程如何随机,但是其期望依然会收敛到一个具体的值上,就像有一只魔法之手在牵引着每次的取值。

下面代码计算𝑛 = 1到 sample_size 时样本的均值:

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt

figsize(12.5, 5)
import pymc as pm

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)
print N_samples

for k in range(3):

    samples = poi(lambda_, size=sample_size)

    partial_average = [samples[:i].mean() for i in N_samples]   # 从[1:sample_size]不同样本数下,三个泊松分布的均值

    plt.plot(N_samples, partial_average, lw=1.5, label="average \
of  $n$ samples; seq. %d" % k)


plt.plot(N_samples, expected_value * np.ones_like(partial_average),
         ls="--", label="true expected value", c="k")

plt.ylim(4.35, 4.65)
plt.title("Convergence of the average of \n random variables to its \
expected value")
plt.ylabel("average of $n$ samples")
plt.xlabel("# of samples, $n$")
plt.legend()

plt.show()

从上图中,我们可以清晰地看出,当 n 很小时,样本均值的变动很大,随着 n 的增大都在逐渐趋向平缓,最终在 4.5 水平线附近小幅摆动。数学家和统计学家将这种“小幅摆动”称为收敛。

我们所熟知的定理:泊松分布的均值等于其参数𝜆。这里所谓的泊松分布指的是一个真实概率分布,但是真实概率分布只有上帝知道我们并不知道,但是当样本数足够大时,根据大数定理,我们可以用一定数量的样本均值来近似代替真实分布的均值,进而得到真实分布的期望。

1. 收敛到期望值的速度

明白了分布的收敛特性之后,另一个与此相关的问题是:收敛到期望值的速度有多快?

1)通过实验方式来观察收敛速度

我们选择一个特定的 n,然后重复实验数千次,并计算与实际期望值的距离的均值:

以上公式可以理解为:前 n 个采样的样本均值与真实值的距离(平均意义上的)

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)

sample_size = 100000
expected_value = lambda_ = 4.5
poi = pm.rpoisson
N_samples = range(1, sample_size, 100)

N_Y = 250  # use this many to approximate D(N)
N_array = np.arange(1000, 50000, 2500)  # use this many samples in the approx. to the variance.
D_N_results = np.zeros(len(N_array))

def D_N(n):
    """
    This function approx. D_n, the average variance of using n samples.
    """
    Z = poi(lambda_, size=(n, N_Y))
    average_Z = Z.mean(axis=0)
    return np.sqrt(((average_Z - expected_value) ** 2).mean())


for i, n in enumerate(N_array):
    D_N_results[i] = D_N(n)


plt.xlabel("$N$")
plt.ylabel("expected squared-distance from true value")
plt.plot(N_array, D_N_results, lw=3,
         label="expected distance between\n\
expected value and \naverage of $N$ random variables.")
plt.plot(N_array, np.sqrt(expected_value) / np.sqrt(N_array), lw=2, ls="--",
         label=r"$\frac{\sqrt{\lambda}}{\sqrt{N}}$")
plt.legend()
plt.title("How 'fast' is the sample average converging? ")

plt.show()

如我们预期,随着 N 增大,样本均值与实际值间距的期望逐渐减小。但也注意到,这个收敛速率也在降低。

样本点从 10000 增加到 20000 时,均值和期望的距离从 0.020 降低到 0.015,减少了 0.005;当样本点从 20000 增加到 40000 时,样本均值和期望之间的距离从 0.015 降低到 0.010,仅仅减少了 0.005。 说明收敛的速率在逐渐变小。

2)通过公式方式来描述收敛速率

实际上收敛的速率是可以通过公式衡量的,在上图的橙色曲线中,是通过下面的公式生成的:

该公式模拟了样本的均值和真实期望之间的距离。

这个特性是非常有趣的,因为对于一个给定的𝑁,就可以知道估计出的样本平均值和真实的期望值之间的误差大小。也很容易直观地想象,如果当前方差很小了,说明摆动幅度已经很小了,基本上肯定是快要进入收敛状态了,所以收敛速率自然也减小了。

0x3:连接概率和频率之间的桥梁 - 大数定理

在开始这个小节讨论之前,我们先抛出一个公式:

频率统计 <---->(大数定理) <----> 期望 <----> 概率

在期望值和估计的概率之间存在一种隐性关系。

首先定义如下函数:

,这种0-1分布有非常好的特性,可以很清晰地表达出期望和概率的关系。

根据大数定理,如果有样本𝑋𝑖 ,可以估计出事件 A 的概率,用𝑃(𝐴)表示

由上面只是函数的定义可知,只有事件发生时1𝐴(𝑥)才取 1 其它情况都取 0, 所以只要用试验中事件发生的总次数除以试验的次数就是概率(回想一下平时我们是怎 么样用频率近似概率的)

在这里,大数定理充当了一个桥梁的作用,将频率统计和概率连接了起来。而连接的水泥就是期望

同时,我们也要明白,大数定理在 N 很小的时候会失效,在 N 很小的时候,此时得到的近似结果是不可靠的。

 

2. 小数据的无序性 - 当 N 很小时大数定理的表现如何?

只有在 N 足够大时,大数定理才成立,然后数据量并非总是足够大的。如果任意运用这一定理,就会导致很荒谬的错误,我们这章来讨论这个问题。

0x1:实例1 - 地理数据的聚合

有些数据会以某一属性进行分组。例如有些数据按照国家,县区或者城市进行分组。当然人口数量也会根据地理区域的不同而不同。如果要对某一地理区域中的一些特征进行统计平均时,必须要注意大数定理适用范围,并了解在人口数较少时它是如何失效的。
假设在数据集 中包含有 5000 个县的数据。除此之外人口数量在整个国家的分布符合 100 到 1500 的均匀分布。 我们关心的是每个县中人口的平均高度。

我们假设不论生活在哪里,每个人的身高都服从相同的分布:height~Normal(150,15)

把数据按照县级别进行分组,所以只能对每个县里的数据做平均。数据如下图:

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 4)
std_height = 15
mean_height = 150

n_counties = 5000
pop_generator = pm.rdiscrete_uniform
norm = pm.rnormal

# generate some artificial population numbers
population = pop_generator(100, 1500, size=n_counties)

average_across_county = np.zeros(n_counties)
for i in range(n_counties): # 5000个县
    # generate some individuals and take the mean
    average_across_county[i] = norm(mean_height, 1. / std_height ** 2,  # 每个县都服从相同的分布: height~Normal(150,15)
                                    size=population[i]).mean()  # 不同的县人口

# located the counties with the apparently most extreme average heights.
i_min = np.argmin(average_across_county)
i_max = np.argmax(average_across_county)

# plot population size vs. recorded average
plt.scatter(population, average_across_county, alpha=0.5, c="#7A68A6")
plt.scatter([population[i_min], population[i_max]],
            [average_across_county[i_min], average_across_county[i_max]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="extreme heights")

plt.xlim(100, 1500)
plt.title("Average height vs. County Population")
plt.xlabel("County Population")
plt.ylabel("Average height in county")
plt.plot([100, 1500], [150, 150], color="k", label="true expected \
height", ls="--")
plt.legend(scatterpoints=1)

plt.show()

从上图中可以观察到什么现象?如果不考虑人口数量大小,推断出的结果就可能存在错误。

例如,如果不考虑人口数大小,我们可能会说有最高身高人口的县也有最低身高的人口,反之亦然(不考虑人口数量,图中的 x 轴就变成一个点了,红色圆圈的点就会在一条竖线上了,所以就反映出最高人口的县也会有最低人口,最低人口的县也会有最高人口)。

1. 在小数据集情况下大数定理失效了吗?

但是这种推断结果真的是对的吗?

从贝叶斯不确定性理论的角度出发,准确的答案应该是:不确定!

理论上来说,确实存在可能在说图中红圈所在的 x 轴对应的县确实存在极值的身高,但是在这种小数据情况下,不确定性是非常高的,高到我们几乎不能将其采纳作为一种推断结果。换句话说,这种高风险的推断对我们来说是没有意义的。

从大数定理的角度来说,这个县不一定有最极端高度的人存在。错误的原因是人口数量少时计算的均值并不能反映出人口真实的期望值(真实是 μ = 150)。样本的数量 N 越小,大数定理就会失效的越厉害。

0x2:实例2 - Kaggle的美国人口普查反馈比例预测

下面是 2010 年美国人口普查数据,它不是按照县对数据分区,而是按照街区进行分区。

本例中的数据来自于 Kaggle 机器学习竞赛中使用的数据集。竞赛的目的是为了预测各个街区对人口普查邮件的回复率,回复率的取值范围为 0 到 100,使用的预测特征包括:人口平均收入,女性数量,停车场的个数,小孩的平均个数等。

下面画出人口调查邮件的回复率街区中人口数量的关系图(多维自变量不方便在一张图上展示)

# -*- coding: utf-8 -*-

import numpy as np
from IPython.core.pylabtools import figsize
import matplotlib.pyplot as plt
import pymc as pm

figsize(12.5, 6.5)
data = np.genfromtxt("./data/census_data.csv", skip_header=1,
                     delimiter=",")
plt.scatter(data[:, 1], data[:, 0], alpha=0.5, c="#7A68A6")
plt.title("Census mail-back rate vs Population")
plt.ylabel("Mail-back rate")
plt.xlabel("population of block-group")
plt.xlim(-100, 15e3)
plt.ylim(-5, 105)

i_min = np.argmin(data[:, 0])
i_max = np.argmax(data[:, 0])

plt.scatter([data[i_min, 1], data[i_max, 1]],
            [data[i_min, 0], data[i_max, 0]],
            s=60, marker="o", facecolors="none",
            edgecolors="#A60628", linewidths=1.5,
            label="most extreme points")

plt.legend(scatterpoints=1)

plt.show()

上图反映了典型的统计现象。这里的典型指的是散点图的形状。散点图是一个典型的三 角形,随着样本数量的增加,这个三角形变得越紧实(因为越来越符合大数定理的要求)。

可以这么理解,随着 N 不断增大,后验概率越来越收敛到某一个固定值,而在 N 很小的时候,后验概率的摆动是很剧烈的。

0x3:实例3 - 如何对Reddit网站上的评论进行排序

上面的两个例子都旗帜鲜明地表达了一个观点:大数据!big data!概率统计和机器学习只有在大数据时才能体现出威力,如果你没有大数据,那就别痴心妄想做机器学习项目了。

但是世上之事,不如意者十之,我们不可能永远在项目中能拿到大数据,就算你身在BAT、google这样的大数据公司,你所拥有的也不一定就是大数据。

参考Andrew Gelman的一句名言:样本从来都不是足够大的。如果 N 太大不足以进行足够精确的估计,你需要获得更多的数据。但当 N “足够大”,你可以开始通过划分数据研究更多的问题,例如在民意调查中,当你已经对全国的民意有了较好的估计,你可以开始分性别、地域、年龄进行更多的统计。N 从来都无法做到足够大,因为当它一旦大了,你总是可以开始研究下一个问题从而需要更多的数据。

回到我们这个小节的问题考虑在线商品的评分问题,这在生活中很常见,各种商城类app里都会有针对某个商品的历史购买者评论。

你如何评价只有一个人,两个人或者三个人对一个商品的五星级评价结果。我们多少都会明白,这样少的评分结果无法真正的反映出商品的真实价值。这种评价信息不充分的情况给排序带来了困难。很多人都知道在线搜索图书, 视频时,如果用它们的评分或者评论对返回结果进行排序,最终用户得到的结果并不会很好。通常情况下搜索结果中排在最前面的视频或评论往往有很好的评分,而这些好的评分都是由粉丝评出的。还有一些真正高质量的视频或者评论排在了稍微靠后的位置, 它们的评分在 4.8 左右。我们如何对这种情况进行修正?

1. 我们如何对评论进行排序

Reddit是一个流行的网站。它主要提供一些小说或者图片的链接,这个网站最流行的部分是它对每个链接的评论。Reddit可以对每个评论进行投票(被称为赞成票或否决票)。默认状态下Reddit将把最好的评论排在前面。

现在假设我们是Reddit后台的算法工程师,我们会如何设计这种排序算法呢?

需要明白的是这是一个非常复杂且牵一发而动全身的事情,一个简单的例子是消费者公众是具有从众效应的,即大多数时候我们会优先看大多数的意见,且很容易接受所谓的大多数人的意见。这种现象表现在很多方面,例如百度搜索的首页的访问量会远大于后面的页面,且这种差距一旦建立,会越来越大;一个商品的评论一旦上了淘宝的爆款热搜,它的购买热度会越来越大,因为公众会优先购买它们最先看到的所谓爆款。

而评论排序这个问题也是一样的,一旦一个评论被排到靠前的位置,用户会优先看到这些评论,而有更大的几率接受该评论的观点,从而形成从众效应,最终导致靠前的评论的点赞数会越来越多。

基于这样的复杂性原因,我们接下来对评论排序这件事展开多维度的讨论。

1)热度

认为点赞数越多的评论越好。但是一个问题在于,当一评论有几百个赞,但是被反对了上千次时,虽然热度很高,但是把这条评论当成最佳评论是有争议的。

但是这里的问题也非常复杂,一个观点有上百文赞同,有上千人反对,那么这个观点就一定是对,或者一定就是错的吗?都不一定,很多时候,这取决于下游的商业和市场策略。

2)差数

用点赞数减去反对数来打分。这解决了以上热度指标的问题。但是依然有缺陷,它不能解决评论的时间属性导致的问题。

因为一个链接在线上会存在很多,不断有用户来点赞或反对,更老的评论会随着时间积累更多的点赞数,但并不意味着他们是最佳评论。这个问题的根源显然是没有进行归一化。

但是反过来想,这个问题也十分复杂,为什么一定要归一化呢?试想,一个评论如果在一个相当长的时间区间内,持续不断地被其他用户所接受并点赞,这不就意味着这个评论的观点接受度很高吗?

3)时间调节

考虑用差数除以评论的寿命(相当于归一化)。这样得到了一种类似每秒差数或每分钟差数这样的速率值。

看起来足够合理了,但是同样存在缺陷,比如一个一秒前发布的评论,只需要一票,就能击败一百秒前的有用九十九票的评论。出现这种问题的原因是没有把持续时间的跨度考虑进去。

避免这一问题的一种解法是只针对 t 秒前的评论进行统计(过滤掉刚刚出现的新鲜的评论)。但是 t 怎么选呢?而且这样是不是意味着 t 秒内的评论都是不好的呢?

读者看出来了吗?我们是在对不稳定量和稳定量进行比较(新评论 PK 老评论)

4)好评率

对某一评论的好评数除以对该评论的好评和差评的总和,得到好评率, 作为排名权重。这解决了 3 中的问题,如果用好评率作为排名结果,即使一个新的评论也会和老的评论一样得到较高的排名得分。

但是依然存在缺陷,在这种度量标准下,如果对一个评论的好评个数只有一个,即好评率为 1.0,而另一个评论的好评个数为 999,差评个数为 1,此时这个评论的好评率为 0.999,显然后一个评论的质量要比前一个高,这“可能”不太客观。

我说“可能”是有理由的,前者虽然只有一个赞成票,但还是有可能由于有 999 个赞成票的后者,只是说这个可能性有多少罢了,我们不能武断地下结论。

2. 真实的好评率如何得到? - 贝叶斯推断来帮忙

我们真正目标就是估计真实的好评率。注意,真实的好评率并不是观测到的好频率, 真实的好评率隐藏起来了,我们能观测到的仅仅是每个评论的好评和差评。

当有 999 个 好评,1 个差评时,我们确实可以说真实的好评率接近于 1,我们能这样说的原因是因为有大数定理的存在;

反之,当只有一次好评时,我们无法断言真实的好评率就是 1,因为在小数据下,后验估计时不确定的,这听起来像个贝叶斯问题,我们显然无法得到一个对好评率的准确估计,我们只能对好评率可能的概率分布进行一个贝叶斯推断。

贝叶斯方法需要好评率的先验知识,获得先验知识的一个方法是基于好评率的历史分布。可以用好评率的历史分布作为定好评率的先验分布。只要对Reddit的评论进行挖掘 就可以得到好评率的历史分布。计算好评率的历史分布时还应注意以下问题:

1. 倾斜数据: 大部分的评论都没有投票数据,因此有许多评论的好评率都非常小 (极端值),这就使分布大多数集中在极端值附件。可以为投票的个数设个阈值,大于 该阈值的情况才予以考虑。阈值的选择还要考虑其对好评率精度的影响。
2. 偏差数据: Reddit有许多不同的子页面,称为subreddits。例如有两个页面,第一个为:r/aww 页面粘贴许多可爱动物的照片,第二个为:r/politics 讨论政治问题。用户对这两个页面的评论行为是不 相同的。浏览第一个页面的人的性格会更友好,更有爱心,也就会有更多的好评(与第二个页面相比),而第二个页面的评论可能会有更多的争论或不同的观点出现。因此并不是所有的评论都是一致的。

鉴于此,我觉得用均匀分布作为评论好评率的先验分布是最佳的做法。

假设共有𝑁个投票数,好评率为𝑝,好评的个数就是二项随机变量(好评或差评),参数为 𝑝 和𝑁。 构造一个函数对 𝑝 进行贝叶斯推断。

import pymc as pm


def posterior_upvote_ratio(upvotes, downvotes, samples=20000):
    """
    This function accepts the number of upvotes and downvotes a particular submission received, 
    and the number of posterior samples to return to the user. Assumes a uniform prior.
    """
    N = upvotes + downvotes
    upvote_ratio = pm.Uniform("upvote_ratio", 0, 1)
    observations = pm.Binomial("obs", N, upvote_ratio, value=upvotes, observed=True)
    # do the fitting; first do a MAP as it is cheap and useful.
    map_ = pm.MAP([upvote_ratio, observations]).fit()
    mcmc = pm.MCMC([upvote_ratio, observations])
    mcmc.sample(samples, samples / 4)
    return mcmc.trace("upvote_ratio")[:]

下面是后验分布的结果:

figsize(11., 8)
posteriors = []
colours = ["#348ABD", "#A60628", "#7A68A6", "#467821", "#CF4457"]
for i in range(len(submissions)):
    j = submissions[i]
    posteriors.append(posterior_upvote_ratio(votes[j, 0], votes[j, 1]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.9,
             histtype="step", color=colours[i % 5], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=18, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )

plt.legend(loc="upper left")
plt.xlim(0, 1)
plt.title("Posterior distributions of upvote ratios on different submissions")

从上图中可以看出,某些分布很窄,其他分布表现为长尾。这体现了我们对真实好评率的不确定性。

3. 如何针对评论的后验分布的进行排序

回到我们这个小节的主题,我们的目标是:如何对评论进行从好到坏的排序?

当然,我们是无法对分布进行排序的,排序只能是标量值。有很多种方法能够从分布中提取出标量值,用期望或均值来表示分布就是一个方法,但是均值并不是一个好办法,因为它没有考虑到分布的不确定性。

我们这里建议使用 95%最小可信值,定义为真实参数只有 5% 的可能性低于该值(即 95% 置信区间的下界),下图是根据 95%最小可信值画的后验分布。

N = posteriors[0].shape[0]
lower_limits = []

for i in range(len(submissions)):
    j = submissions[i]
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.9,
             histtype="step", color=colours[i], lw=3,
             label='(%d up:%d down)\n%s...' % (votes[j, 0], votes[j, 1], contents[j][:50]))
    plt.hist(posteriors[i], bins=20, normed=True, alpha=.2,
             histtype="stepfilled", color=colours[i], lw=3, )
    v = np.sort(posteriors[i])[int(0.05 * N)]
    # plt.vlines( v, 0, 15 , color = "k", alpha = 1, linewidths=3 )
    plt.vlines(v, 0, 10, color=colours[i], linestyles="--", linewidths=3)
    lower_limits.append(v)
    plt.legend(loc="upper left")

plt.legend(loc="upper left")
plt.title("Posterior distributions of upvote ratios on different submissions");
order = np.argsort(-np.array(lower_limits))
print(order, lower_limits)

为何基于这个量进行排序是个好主意?

因为基于 95% 最小可信值进行排序,是对最佳评论的一个最为保守的估计。即便在最差情况下,也就是说我们明显高估了好评率时,也能确保最佳的评论排在顶端。在这一排序思路中,我们利用了以下很自然的特性:

1. 给定两个好评率相同的评论,我们会选择票数最多的作为更加评论,因为票数更多的评论在贝叶斯后验概率分布中更集中,落在 95%置信区间的概率也更大,也即我们更确定其好评率更好。
2. 给定两个票数一样的评论,我们会选择好评数更多的。

 

3. 大数定理的一些更多思考

尽管大数定理非常酷,但是正如其名,它只适用于足够大的数据量。我们已经看到,如果不考虑数据的构造,那么估计结果会受到很大影响。

1. 通过(低成本地)获得大量后验样本,可以确保大数定理适用于期望的近似估计
2. 在样本数量很小的贝叶斯推断中,后验分布包括更多的随机性。从后验分布中将能看出这些随机性的存在。后验分布越平坦随机性越大,越陡峭随机性越小。因此推断的结果是需要调校的。
3. 不考虑样本量会带来很大的影响,此时排序的依据往往是不稳定的,并会导致非常病态的排序结果。

 

posted @ 2018-09-02 11:27 LittleHann 阅读(...) 评论(...) 编辑 收藏