第3章 概率统计基础 概率统计知识在人工智能领域发挥着非常重要的作用，如深度学习理论、概率图模型等都依赖于概率分布作为框架的基本建模语言，与此同时，复杂的深度学习等算法又超越了传统统计分析方法，在人工智能与大数据领域发挥了不可替代的作用。本章将从概率论与数理统计的基本概念和原理出发，重点介绍与人工智能有关的数理统计方法，并以Python语言为例设计对应的实验。 3．1概率论 3．1．1概率与条件概率 在我们的日常生活和工作中，经常会遇到许多包含不确定结果的现象，如股票价格是涨还是跌呢？已经布满乌云的天空是否会下雨呢？某项工程是否能够按期完成呢？投资项目盈利可能性有多大呢？上述现象都有两个特点： ①事先不能确定哪一个结果会出现； ②各种结果在多次重复过程中可能会体现某种规律。 人们把这类现象称为随机现象或者不确定现象，同时使用一个数值来度量随机现象中某一结果出现可能性的大小，这个数值就被称为概率(Probability)。概率的取值在0到1之间，概率值等于0表明该现象不可能发生，概率值等于1表明该现象必然发生，介于0和1之间的概率则说明该现象出现可能性的不同程度。 随机试验是对随机现象的观察、记录、实验的统称，是在相同条件下对某随机现象进行的大量重复观测，具有以下特点： ① 在试验前不能断定其将发生什么结果，但可明确指出或说明试验的全部可能结果是什么； ② 在相同的条件下试验可大量地重复； ③ 重复试验的结果是以随机方式或偶然方式出现。 人们把随机试验E的所有可能结果组成的集合称为E的样本空间，记为S。样本空间的元素即E的每个结果称为样本点。 假设在相同条件下，进行了n次试验，在这n次试验中，人们把事件A发生的次数称为事件A发生的频数，A发生的次数与试验次数的比值称为事件A发生的频率。 当试验次数增加时，随机事件A发生的频率趋于一个稳定值，记为p，p就称为该事件发生的概率，记为P(A)=p。 条件概率(Conditional Probability)是一种带有附加条件的概率，例如，如果事件A与事件B是相依事件，即事件A的概率随事件B是否发生而变化，记为P(A|B)，表示在事件B发生的条件下，事件A发生的概率，相当于A在B中所占的比例。 【例31】一个家庭中有两个小孩，已知至少一个是女孩，问两个都是女孩的概率是多少(假定生男生女的可能性是相等的)？ 解： 由题意，样本空间为 S={(兄，弟)，(兄，妹)，(姐，弟)，(姐，妹)} B={(兄，妹)，(姐，弟)，(姐，妹)} A={(姐，妹)} 由于事件B已经发生，所以这时试验的所有可能只有3种，而事件A包含的基本事件只占其中的一种，所以有P(A|B)=1/3，即在已知至少一个是女孩的情况下，两个都是女孩的概率为1/3。 在这个例子中，如果不知道事件B发生，则事件A发生的概率为P(A)=1/4。这里P(A)≠P(A|B)，其原因在于事件B的发生改变了样本空间，使它由原来的S缩减为新的样本空间SB=B。 3．1．2随机变量 随机变量(Random Variable)表示随机试验各种结果的实值单值函数，即能用数学分析方法来研究随机现象。例如某一时间内公共汽车站等车的乘客人数、淘宝在一定时间内的交易次数、百度在一段时间内某个关键词搜索次数等，都是随机变量的实例。这些例子中所提到的量，尽管它们的具体内容各式各样，但从数学观点来看，它们表现了同一种情况，就是每个变量都可以随机地取得不同的数值，而在进行试验或测量之前，我们要预言这个变量将取得某个确定的数值是不可能的。 按照随机变量可能取得的值，可以把它们分为两种基本类型： 离散型和连续型。 离散型随机变量即在一定区间内变量取值为有限个或可数个。例如某地区某年人口的出生数、死亡数，某药治疗某病病人的有效数、无效数等。离散型随机变量根据不同的概率分布有伯努利分布、二项分布、几何分布、泊松分布、超几何分布等。 连续型随机变量即在一定区间内变量取值有无限个，或数值无法一一列举出来。例如某地区男性健康成人的身高值、体重值，一批传染性肝炎患者的血清转氨酶测定值等。连续型随机变量根据不同的概率分布有均匀分布、指数分布、正态分布、伽马分布等。 随机变量的性质主要有两类： 一类是大而全的性质，这类性质可以详细描述所有可能取值的概率，例如描述连续型随机变量的累积分布函数 (Cumulative Distribution Function，CDF)、概率密度函数(Probability Density Function，PDF)，描述离散型随机变量的概率质量分布函数(Probability Mass Function，PMF)等； 另一类是找到该随机变量的一些特征或代表值，例如随机变量的方差(Variance)、期望(Expectation)、置信区间等数字特征。 SciPy、NumPy、Matplotlib是Python中使用最为广泛的科学计算工具包，基本上可以处理大部分的统计与可视化作图等任务。SciPy有个stats模块，这个模块中包含了概率论及统计相关的函数。下面我们用Python来模拟实现各种随机变量的分布及其图形化显示。 3．1．3离散随机变量分布Python实验 1． 伯努利分布 伯努利分布(Bernoulli Distribution)又称两点分布或01分布，其样本空间中只有两个点，一般取为{0，1}，不同的伯努利分布只是取到这两个值的概率不同。伯努利分布只有一个参数p，记作X~Bernoulli(p)，或X~B(1，p)，读作X服从参数为p的伯努利分布。 下面以抛硬币为例模拟伯努利分布。如果将抛一次硬币看作一次伯努利实验，且将正面朝上记为1，反面朝上记为0。那么伯努利分布中的参数p就表示硬币正面朝上的概率。 伯努利分布的概率分布用Python代码绘制如下： #第3章/bernoulli_pmf.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def bernoulli_pmf(p=0.0): ber_dist = stats.bernoulli(p) x = ［0, 1］ x_name = ［'0', '1'］ pmf = ［ber_dist.pmf(x［0］), ber_dist.pmf(x［1］)］ plt.bar(x, pmf, width=0.15) plt.xticks(x, x_name) plt.ylabel('Probability') plt.title('PMF of bernoulli distribution') plt.show() bernoulli_pmf(p=0.3) 图31是该程序的运行结果图。Probability为概率，横坐标表示随机变量X，1表示正面朝上。 图31伯努利分布B(1，0．3)的PMF柱状图 通常为了得到比较准确的某个服从伯努利分布的随机变量的期望，需要大量重复伯努利试验，例如重复n次，然后利用“正面朝上的次数/n”来估计p。 2． 二项分布 如果把一个伯努利分布独立地重复n次，就得到了一个二项分布。二项分布是最重要的离散型概率分布之一。随机变量X要满足这个分布有两个重要条件： ① 各次试验的条件是稳定的； ② 各次试验之间是相互独立的。 还是以利用抛硬币的例子来比较伯努利分布和二项分布的区别。如果将抛一次硬币看作一次伯努利实验，且将正面朝上记为1，反面朝上记为0。那么抛n次硬币，记录正面朝上的次数Y，Y就服从二项分布。假如硬币是均匀的，Y的取值应该大部分集中在n/2附近，而非常大或非常小的值都很少。由此可见，二项分布关注的是计数，而伯努利分布关注的是比值(正面朝上的计数/n)。一个随机变量X服从参数为n和p的二项分布，记作X~Binomial(n，p)或X~B(n，p)。 下面利用Python代码模拟抛一枚不均匀的硬币20次，设正面朝上的概率为0．6。 #第3章/binom_dis.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def binom_dis(n=1, p=0.1): binom_dis = stats.binom(n, p) x = np.arange(binom_dis.ppf(0.0001), binom_dis.ppf(0.9999)) print(x)#［ 0.1.2.3.4.］ fig, ax = plt.subplots(1, 1) ax.plot(x, binom_dis.pmf(x), 'bo', label='binom pmf') ax.vlines(x, 0, binom_dis.pmf(x), colors='b', lw=5, alpha=0.5) ax.legend(loc='best', frameon=False) plt.ylabel('Probability') plt.title('PMF of binomial distribution(n={}, p={})'.format(n, p)) plt.show() binom_dis(n=20, p=0.6) 以上定义了一个n=20，p=0．6的二项分布，表示每次试验抛硬币，该硬币正面朝上的概率大于背面朝上的概率，共抛20次并记录正面朝上的次数。该分布的概率质量分布函数如图32所示，Probability为概率，X表示随机变量： 正面朝上的次数。 图32二项分布B(20，0．6)的PMF图 从图32中可以看出，该分布的概率质量分布函数PMF图明显向右边偏移，在x=12处取到最大概率，这是因为这个硬币正面朝上的概率大于反面朝上的概率。 为了比较准确地得到某个服从二项分布的随机变量的期望，需要大量重复二项分布试验，例如有m个人进行试验，每人抛n次，然后利用“所有人得到的正面次数之和/m”来估计n和p，总共相当于做了n×m次伯努利实验。 3． 泊松分布 日常生活中很多事件的发生是有固定频率的，例如某医院平均每小时出生的婴儿数，某网站平均每分钟的访问数、某一服务设施在一定时间内收到的服务请求的次数、汽车站台的候客人数等。它们的特点是可以预估这些事件在某个时间段内发生的总次数，但是没法知道具体的发生时间。如果某事件以固定强度λ随机且独立地出现，该事件在单位时间内出现的次数(个数)可以看成是服从泊松分布。我们把一个随机变量X服从参数为λ的泊松分布，记作X~Poisson(λ)，或X~P(λ)。泊松分布适合于描述单位时间内随机事件发生次数的概率分布。 下面是参数μ=8时的泊松分布Python实现，在SciPy中将泊松分布的参数表示为μ。 #第3章/poisson_pmf.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def poisson_pmf(mu=3): poisson_dis = stats.poisson(mu) x = np.arange(poisson_dis.ppf(0.001), poisson_dis.ppf(0.999)) print(x) fig, ax = plt.subplots(1, 1) ax.plot(x, poisson_dis.pmf(x), 'bo', ms=8, label='poisson pmf') ax.vlines(x, 0, poisson_dis.pmf(x), colors='b', lw=5, alpha=0.5) ax.legend(loc='best', frameon=False) plt.ylabel('Probability') plt.title('PMF of poisson distribution(mu={})'.format(mu)) plt.show() poisson_pmf(mu=8) 代码运行后输出的概率质量分布PMF如图33所示。 图33二项分布B(20，0．6)的PMF图 如果仅仅看二项分布与泊松分布的概率质量分布图，也可以发现它们的相似度非常高。泊松分布可以作为二项分布的极限。一般来说，若X~B(n，p)，其中n很大，p很小，而np=λ不太大时，则X的分布接近于泊松分布P(λ)。 3．1．4连续随机变量分布Python实验 1． 均匀分布 均匀分布(Uniform Distribution)是最简单的连续型概率分布，因为其概率密度是一个常数，不随随机变量X取值的变化而变化。如果连续型随机变量 X 具有如下的概率密度函数，则称X 服从 ［a，b］ 上的均匀分布，记作 X~U(a，b)或 X~Unif(a，b)。 fX(x)=1b－aab(31) 从式(31)可以看出，定义一个均匀分布需要两个参数，定义域区间的起点a 和终点 b，在Python中用 location 和 scale分别表示起点和区间长度，代码如下： #第3章/uniform_dis.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def uniform_distribution(loc=0, scale=1): uniform_dis = stats.uniform(loc=loc, scale=scale) x = np.linspace(uniform_dis.ppf(0.01), uniform_dis.ppf(0.99), 100) fig, ax = plt.subplots(1, 1) #直接传入参数 ax.plot(x, stats.uniform.pdf(x, loc=2, scale=4), 'r-', lw=5, alpha=0.6, label='uniform pdf') #从冻结的均匀分布取值 ax.plot(x, uniform_dis.pdf(x), 'k-', lw=2, label='frozen pdf') #计算ppf分别等于0.001、0.5和0.999时的x值 vals = uniform_dis.ppf(［0.001, 0.5, 0.999］) print(vals)#［ 2.0044.5.996］ #检测cdf 和 ppf的精确度 print(np.allclose(［0.001, 0.5, 0.999］, uniform_dis.cdf(vals)))#结果为True r = uniform_dis.rvs(size=10000) ax.hist(r, normed=True, histtype='stepfilled', alpha=0.2) plt.ylabel('Probability') plt.title(r'PDF of Unif({}, {})'.format(loc, loc+scale)) ax.legend(loc='best', frameon=False) plt.show() uniform_distribution(loc=2, scale=4) 上面的代码采用两种方式： 直接传入参数和先冻结一个分布，画出来均匀分布的概率分布函数，此外还从该分布中取了10000个值作直方图，运行结果如图34所示。 图34均匀分布 U(2，6)的PDF图 2． 指数分布 指数分布和离散型的泊松分布之间有很大的关系。泊松分布表示单位时间(或单位面积)内随机事件的平均发生次数，指数分布则可以用来表示独立随机事件发生的时间间隔。由于发生次数只能是自然数，所以泊松分布自然就是离散型的随机变量，而时间间隔则可以是任意的实数，因此其定义域是(0，+∞)。 如果一个随机变量X的概率密度函数满足以下形式，就称X为服从参数λ的指数分布(Exponential Distribution)，记作X~E(λ)或 X~Exp(λ)。指数分布只有一个参数λ，且λ>0。 fX(x)=λe－λex>0 0其他(32) 指数分布的一个显著的特点是其具有无记忆性。例如，如果排队的顾客接受服务的时间长短服从指数分布，那么无论你已经排了多长时间的队，在排t分钟的概率始终是相同的。用公式表示为P(X≥s+t|X≥s)=P(X≥t)，代码如下： #第3章/exponential_dis.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def exponential_dis(loc=0, scale=1.0): """ 指数分布，exponential continuous random variable 按照定义，指数分布只有一个参数lambda，这里的scale = 1/lambda :param loc:定义域的左端点，相当于将整体分布沿x轴平移loc :param scale: lambda的倒数，loc+scale表示该分布的均值，scale^2表示该分布的方差 :return: """ exp_dis = stats.expon(loc=loc, scale=scale) x = np.linspace(exp_dis.ppf(0.000001), exp_dis.ppf(0.999999), 100) fig, ax = plt.subplots(1, 1) #直接传入参数 ax.plot(x, stats.expon.pdf(x, loc=loc, scale=scale), 'r-', lw=5, alpha=0.6, label='uniform pdf') #从冻结的均匀分布取值 ax.plot(x, exp_dis.pdf(x), 'k-', lw=2, label='frozen pdf') #计算ppf分别等于0.001、0.5和0.999时的x值 vals = exp_dis.ppf(［0.001, 0.5, 0.999］) print(vals)#［ 2.0044.5.996］ #检测cdf 和 ppf的精确度 print(np.allclose(［0.001, 0.5, 0.999］, exp_dis.cdf(vals))) r = exp_dis.rvs(size=10000) ax.hist(r, normed=True, histtype='stepfilled', alpha=0.2) plt.ylabel('Probability') plt.title(r'PDF of Exp(0.5)') ax.legend(loc='best', frameon=False) plt.show() exponential_dis(loc=0, scale=2) Exp(0．5)的PDF如图35所示。 图35二项分布B(20，0．6)的PDF图 3． 正态分布 正态分布(Normal Distribution)，也称常态分布，又名高斯分布(Gaussian Distribution)，最早由A．棣莫弗在求二项分布的渐近公式中得到。C．F．高斯在研究测量误差时从另一个角度导出了它。P．S．拉普拉斯和高斯研究了它的性质。正态分布是一个在数学、物理及工程等领域都非常重要的概率分布，在统计学的许多方面有着重大的影响力。正态曲线呈钟形，两头低，中间高，左右对称因其曲线呈钟形，因此人们又经常称其为钟形曲线。 若随机变量X服从一个数学期望为μ、方差为σ2的正态分布，记为N(μ，σ2)。其概率密度函数为正态分布的期望值μ决定了其位置，其标准差σ决定了分布的幅度。当μ=0，σ=1时的正态分布是标准正态分布。其概率密度函数为 f(x)=1σ2πe－(x－μ)22σ2(33) 当μ=0，σ=1时，该正态分布称为标准正态分布。由于标准正态分布在统计学中的重要地位，它的累积分布函数CDF有一个专门的表示符号： Φ。正态分布，Python代码如下： #第3章/normal_dis.py #绘制正态分布概率密度函数 import numpy as np import matplotlib.pyplot as plt import math u = 0#均值μ u01 = -2 sig = math.sqrt(0.2)#标准差δ x = np.linspace(u - 3 * sig, u + 3 * sig, 50) y_sig = np.exp(-(x - u) ** 2 / (2 * sig ** 2)) / (math.sqrt(2 * math.pi) * sig) print(x) print("=" * 20) print(y_sig) plt.plot(x, y_sig, "r-", linewidth=2) plt.grid(True) plt.show() 该程序运行结果如图36所示。 图36标准正态分布图 正态分布中的两个参数含义如下： 当固定σ，改变μ的大小时，f(x)图形的形状不变，只是沿着x轴作平移变换，因此μ被称为位置参数(决定对称轴的位置)； 当固定μ，改变σ的大小时，f(x)图形的对称轴不变，形状改变，σ越小，图形尖峰越陡峭。σ越大，图形越平坦，因此σ被称为尺度参数，决定曲线的分散程度。 3．2数理统计基础 数理统计(Mathematics Statistics)是一门应用性很强的数学分支，它以概率论为基础，根据试验或者观察得到的数据来研究随机现象，对随机现象的概率p做出一些合理的估计和判断。数理统计包含许多内容，但其核心部分是统计推断，它包括两个基本问题，即参数估计和假设检验。 3．2．1总体和样本 统计学中的概念很多，其中有几个概念是经常用到的，包括总体、样本、参数、统计量、变量等。在数理统计中，称研究对象的全体为总体，组成总体的每个基本单元叫个体。从总体X中随机抽取一部分个体X1，X2，…，Xn，称 X1，X2，…，Xn为取自X的容量为n的样本。 例如，为了研究某厂生产的一批元件质量的好坏，规定使用寿命低于1000小时的为次品，则该批元件的全体就为总体，每个元件就是个体。实际上，数理统计学中的总体是指与总体相联系的某个(或某几个)数量指标X取值的全体。例如，该批元件的使用寿命X的取值全体就是研究对象的总体。显然X是随机变量，这时就称X为总体。 为了判断该批元件的次品率，最精确的办法是取出全部元件，对元件的寿命做试验。然而，寿命试验具有破坏性，即使某些试验是非破坏性的，因此只能从总体中抽取一部分，对n个个体进行试验。试验结果可得到一组数值集合{x1，x2，…，xn}，其中每个xi是第i次抽样观察的结果。由于要根据这些观察结果来对总体进行推断，所以对每次抽样需要有一定的要求，要求每次抽取样本必须是随机的、独立的，这样才能较好地反映总体情况。所谓随机是指每个个体被抽到的机会是均等的，这样抽到的个体才具有代表性。若x1，x2，…，xn相互独立，且每个xi与X同分布，则称x1，x2，…，xn为简单随机样本，简称样本。通常把n称为样本容量。 值得注意的是，样本具有两重性，即当在一次具体的抽样后它是一组确定的数值。但在一般叙述中样本也是一组随机变量，因为抽样是随机的。一般地，用X1，X2，…，Xn表示随机样本，它们取到的值记为x1，x2，…，xn称为样本观测值。 样本作为随机变量有一定的概率分布，这个概率分布称为样本分布。显然，样本分布取决于总体的性质和样本的性质。 3．2．2统计量与抽样分布 1． 统计量 数理统计的任务是采集和处理带有随机影响的数据，或者说收集样本并对之进行加工，以此对所研究的问题做出一定的结论，这一过程称为统计推断。在统计推断中，对样本进行加工整理，实际上就是根据样本计算出一些量，使得这些量能够将所研究问题的信息集中起来。这种根据样本计算出的量就是下面将要定义的统计量，因此，统计量是样本的某种函数。 设X1，X2，…，Xn是总体X的一个简单随机样本，T(X1，X2，…，Xn)为一个n元连续函数，且T中不包含任何关于总体的未知参数，则称T(X1，X2，…，Xn)是一个统计量，称统计量的分布为抽样分布。下面列出几个常用的统计量。 设X1，X2，…，Xn是来自总体X的一个样本，x1，x2，…，xn是这一样本的观测值。统计量定义如下。 样本均值 =1n∑ni=1Xi（34） 样本方差 S2=1n－1∑ni=1(Xi－)2（35） 样本标准差 S=S2=1n－1∑ni=1(Xi－)2（36） 样本k阶(原点)矩 Ak=1n∑ni=1Xkik=1，2，…（37） 样本k阶中心矩 Bk=1n∑ni=1(Xi－)kk=2，3，…（38） 2． 三大抽样分布 统计量是样本的函数，它是一个随机变量。统计量的分布称为抽样分布，在使用统计量进行统计推断时常需知道它的分布。当总体的分布函数已知时，抽样分布是确定的，然而要求出统计量的精确分布比较困难。下面介绍几个常用的统计量分布。  卡方(χ2)分布 设X1，X2，…，Xn是来自总体N(0，1)的样本，则统计量χ2=X21+X22+…+X2n所服从的分布称为自由度为n的χ2分布(χ2distribution)，记为χ2~χ2(n)。χ2(n)分布的概率密度函数为 f(y)=12n2Γn2yn2－1e－y2y>0 0其他（39） 不同参数的卡方分布，Python实现代码如下： #第3章/binom_dis.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def diff_chi2_dis(): """ 不同参数下的卡方分布 :return: """ #chi2_dis_0_5 = stats.chi2(df=0.5) chi2_dis_1 = stats.chi2(df=1) chi2_dis_4 = stats.chi2(df=4) chi2_dis_10 = stats.chi2(df=10) chi2_dis_20 = stats.chi2(df=20) #x1 = np.linspace(chi2_dis_0_5.ppf(0.01), chi2_dis_0_5.ppf(0.99), 100) x2 = np.linspace(chi2_dis_1.ppf(0.65), chi2_dis_1.ppf(0.9999999), 100) x3 = np.linspace(chi2_dis_4.ppf(0.000001), chi2_dis_4.ppf(0.999999), 100) x4 = np.linspace(chi2_dis_10.ppf(0.000001), chi2_dis_10.ppf(0.99999), 100) x5 = np.linspace(chi2_dis_20.ppf(0.00000001), chi2_dis_20.ppf(0.9999), 100) fig, ax = plt.subplots(1, 1) #ax.plot(x1, chi2_dis_0_5.pdf(x1), 'b-', lw=2, label=r'df = 0.5') ax.plot(x2, chi2_dis_1.pdf(x2), 'g-', lw=2, label='df = 1') ax.plot(x3, chi2_dis_4.pdf(x3), 'r-', lw=2, label='df = 4') ax.plot(x4, chi2_dis_10.pdf(x4), 'b-', lw=2, label='df = 10') ax.plot(x5, chi2_dis_20.pdf(x5), 'y-', lw=2, label='df = 20') plt.ylabel('Probability') plt.title(r'PDF of $＼chi^2$ Distribution') ax.legend(loc='best', frameon=False) plt.show() diff_chi2_dis() 图37不同参数下卡方分布的PDF图 当自由度df等于1或2时，函数图像都呈单调递减的趋势； 当df大于或等于3时，呈先增后减的趋势。从定义上来看，df的值只能取正整数，如图37所示。 3． t分布 设X~N(0，1)，Y~χ2(n)，并且X和Y独立，则称随机变量t=XYn服从自由度为n的t分布(tdistribution)，记为t~t(n)。 t(n)分布的概率密度函数为 h(t)=Γ(n+1)/2nπΓ(n/2)1+t2n－(n+1)/2－∞0 0其他（311） F分布经常被用来对两个样本方差进行比较。它是方差分析的一个基本分布，也被用于回归分析中的显著性检验。 F分布有两个参数： dfn和dfd，分别代表分子上的第一自由度和分母上的第二自由度。下面是不同参数下F分布的Python实现，代码如下： #第3章/binom_dis.py import numpy as np from scipy import stats import matplotlib.pyplot as plt def diff_f_dis(): """ 不同参数下的F分布 :return: """ #f_dis_0_5 = stats.f(dfn=10, dfd=1) f_dis_1_30 = stats.f(dfn=1, dfd=30) f_dis_30_5 = stats.f(dfn=30, dfd=5) f_dis_30_30 = stats.f(dfn=30, dfd=30) f_dis_30_100 = stats.f(dfn=30, dfd=100) f_dis_100_100 = stats.f(dfn=100, dfd=100) #x1 = np.linspace(f_dis_0_5.ppf(0.01), f_dis_0_5.ppf(0.99), 100) x2 = np.linspace(f_dis_1_30.ppf(0.2), f_dis_1_30.ppf(0.99), 100) x3 = np.linspace(f_dis_30_5.ppf(0.00001), f_dis_30_5.ppf(0.99), 100) x4 = np.linspace(f_dis_30_30.ppf(0.00001), f_dis_30_30.ppf(0.999), 100) x6 = np.linspace(f_dis_30_100.ppf(0.0001), f_dis_30_100.ppf(0.999), 100) x5 = np.linspace(f_dis_100_100.ppf(0.0001), f_dis_100_100.ppf(0.9999), 100) fig, ax = plt.subplots(1, 1, figsize=(20, 10)) #ax.plot(x1, f_dis_0_5.pdf(x1), 'b-', lw=2, label=r'F(0.5, 0.5)') ax.plot(x2, f_dis_1_30.pdf(x2), 'g-', lw=2, label='F(1, 30)') ax.plot(x3, f_dis_30_5.pdf(x3), 'r-', lw=2, label='F(30, 5)') ax.plot(x4, f_dis_30_30.pdf(x4), 'm-', lw=2, label='F(30, 30)') ax.plot(x6, f_dis_30_100.pdf(x6), 'c-', lw=2, label='F(30, 100)') ax.plot(x5, f_dis_100_100.pdf(x5), 'y-', lw=2, label='F(100, 100)') plt.ylabel('Probability') plt.title(r'PDF of f Distribution') ax.legend(loc='best', frameon=False) plt.savefig('f_diff_pdf.png', dip=500) plt.show() diff_f_dis() 不同参数下F分布的PDF图如图39所示。 图39不同参数下F分布的PDF图 3．2．3大数定律与中心极限定理 1． 大数定律 大数定律是一种描述当试验次数很大时所呈现的概率性质的定律，它由概率统计定义“频率收敛于概率”引申而来。简而言之，就是n个独立分布的随机变量其观察值的均值依概率收敛于这些随机变量所属分布的理论均值，也就是总体均值。大数定律为这种后验认识世界的方式提供了坚实的理论基础。统计学中常采用大数定律用样本均值来估计总体的期望。例如，假设每次从1、2、3当中随机选取一个数字，随着抽样次数的增加，样本均值越来越趋近于总体期望((1+2+3)/3=2)。 下面用程序模拟抛硬币的过程来辅助说明大数定律： 用random模块生成区间［0，1)之间的随机数，如果生成的数小于0．5，就记为硬币正面朝上，否则记为硬币反面朝上。由于random．random()生成的数可以看作服从区间［0，1)上的均匀分布，所以以0．5为界限，随机生成的数中大于0．5或小于0．5的概率应该是相同的，相当于硬币出现正反面的概率是均匀的。这样就用随机数模拟出了实际的抛硬币试验。理论上试验次数越多，即抛硬币的次数越多，正反面出现的次数之比越接近于1，也就是说正反面各占一半。 #第3章/binom_dis.py import random import matplotlib.pyplot as plt def flip_plot(minExp, maxExp): """ Assumes minExp and maxExp positive integers; minExp < maxExp Plots results of 2**minExp to 2**maxExp coin flips """ #两个参数的含义，抛硬币的次数为2的minExp次方到2的maxExp次方，也就是一共做了 #(2**maxExp - 2**minExp)批次实验，每批次重复抛硬币2**n次 ratios = ［］ xAxis = ［］ for exp in range(minExp, maxExp + 1): xAxis.append(2**exp) for numFlips in xAxis: numHeads = 0#初始化，硬币正面朝上的计数为0 for n in range(numFlips): if random.random() < 0.5:#random.random()从［0, 1)随机地取出一个数 numHeads += 1#当随机取出的数小于0.5时，正面朝上的计数加1 numTails = numFlips - numHeads#得到本次试验中反面朝上的次数 ratios.append(numHeads/float(numTails))#正反面计数的比值 plt.title('Heads/Tails Ratios') plt.xlabel('Number of Flips') plt.ylabel('Heads/Tails') plt.plot(xAxis, ratios) plt.hlines(1, 0, xAxis［-1］, linestyles='dashed', colors='r') plt.show() flip_plot(4, 16) 该程序顺利执行后将会出现如图310所示界面。 图310大数定律模拟输出 从图310可以看出，随着实验次数的增加，硬币正反面出现次数之比越来越接近于1。 2． 中心极限定理 大数定律揭示了大量随机变量的平均结果，但没有涉及随机变量的分布问题。而中心极限定理说明的是在一定条件下，大量独立随机变量的平均数是以正态分布为极限的。中心极限定理指出大量的独立随机变量之和具有近似于正态的分布。因此，它不仅提供了计算独立随机变量之和的近似概率的简单方法，而且有助于解释为什么有很多自然群体的经验频率呈现出钟形(即正态)曲线这一事实，因此中心极限定理这个结论使正态分布在数理统计中具有很重要的地位，也使正态分布有了广泛的应用。中心极限定理有辛钦中心极限定理、德莫佛拉普拉斯中心极限定理、李亚普洛夫中心极限定理、林德贝尔格定理等表现形式。 下面使用Python来模拟并验证中心极限定理。假设我们现在观测一个人掷骰子。这个骰子是公平的，也就是说掷出1~6的概率都是相同的： 1/6。他掷了一万次，用Python来模拟投掷过程，代码如下： #第3章/central_limit.py import numpy as np random_data = np.random.randint(1, 7, 10000) print random_data.mean()#打印平均值 print random_data.std()#打印标准差 生成的平均值： 3．4927(注意： 每次执行输出结果会略有不同)，生成的标准差： 1．7079。 平均值接近3．5很好理解，因为每次掷出来的结果是1、2、3、4、5、6。 每个结果的概率是1/6，所以加权平均值就是3．5。接下来随机抽取一组抽样，例如从生成的数据中随机抽取10个数字： sample1 = ［］ for i in range(0, 10): sample1.append(random_data［int(np.random.random() * len(random_data))］) print sample1#打印出来 这10个数字的结果是： ［3，4，3，6，1，6，6，3，4，4］。 平均值： 4．0。 标准差： 1．54。 可以看到，只抽10个样本的时候，样本的平均值(4．0)距离总体的平均值(3．5)有所偏差。 下面抽取1000组，每组50个样本，并且把每组的平均值都算出来。 samples = ［］ samples_mean = ［］ samples_std = ［］ for i in range(0, 1000): sample = ［］ for j in range(0, 50): sample.append(random_data［int(np.random.random() * len(random_data))］) sample_np = np.array(sample) samples_mean.append(sample_np.mean()) samples_std.append(sample_np.std()) samples.append(sample_np) samples_mean_np = np.array(samples_mean) samples_std_np = np.array(samples_std) print samples_mean_np 共1000组平均值大概是这样的： ［3．44，3．42，3．22，3．2，2．94 … 4．08，3．74］，结果输出如下： 平均值： 3．48494。 标准差： 0．23506。 然后把这1000组数字用直方图画出来，如图311所示。 图311大数定律下中心极限的正态分布图 从图311中可以看出，其结果已经非常接近正态分布了。 3．3参数估计 参数估计是数理统计研究的主要问题之一。假设总体X~N(μ，σ2)，μ和σ2是未知参数，X1，X2，…，Xn是来自X的样本，样本值是x1，x2，…，xn，要由样本值来确定μ和σ2的估计值，这就是参数估计问题，参数估计分为点估计(Point Estimation)和区间估计(Interval Estimation)。 3．3．1点估计 所谓点估计是指把总体的未知参数估计为某个确定的值或在某个确定的点上，所以点估计又称为定值估计。 设总体X的分布函数为F(x，θ)，θ是未知参数，X1，X2，…，Xn是X的一样本，样本值为x1，x2，…，xn，构造一个统计量(X1，X2，…，Xn)，用它的观察值 (x1，x2，…，xn)作为θ的估计值，这种问题称为点估计问题。习惯上称随机变量(X1，X2，…，Xn)为θ的估计量，称(x1，x2，…，xn)为估计值。 构造估计量(X1，X2，…，Xn)的方法很多，下面介绍常用的矩估计法和极大似然估计法。 1． 矩估计法 矩估计法是一种古老的估计方法。矩是描述随机变量的最简单的数字特征，样本来自于总体，样本矩在一定程度上也反映了总体矩的特征，且在样本容量n增大的条件下，样本的k阶原点矩Ak=1n∑ni=1Xki依概率收敛到总体X的k阶原点矩μk=E(Xk)，即Akpμk(n→∞)，k=1，2，…，因而自然想到用样本矩作为相应总体矩的估计量，而以样本矩的连续函数作为相应总体矩的连续函数的估计量，这种估计方法就称为矩估计法。 矩估计法的一般做法： 设总体X~F(X； θ1，θ2，…，θl)其中θ1，θ2，…，θl均未知。 (1) 如果总体X的k阶矩μk=E(Xk)(1≤k≤l)均存在，则 μk=μk(θ1，θ2，…，θl)，(1≤k≤l) (2) 令μ1(θ1,θ2,…,θl)A1 μ2(θ1,θ2,…,θl)A2  μi(θ1,θ2,…,θl)Al 其中，Ak(1≤k≤l)为样本k阶矩。 求出方程组的解θ^1,θ^2,…,θ^l，称θ^k=θ^k(X1,X2,…,Xn)为参数θk(1≤k≤l)的矩估计量，θ^k=θ^k(x1,x2,…,xn)为参数θk的矩估计值。 【例32】设总体X的密度函数为 f(x； θ)=e－(x－θ)x>θ 0x<θ 其中θ未知，样本为(X1，X2，…，Xn)，求参数θ的矩估计值。 解： A1=。由μ1=A1及 μ1=E(X)=∫+∞－∞xf(x； θ)dx=∫∞0xe－(x－θ)dx=θ+1 有θ=μ1－1，得A1=代替上式中的μ1，得到参数θ的矩估计量为 θ^=－1 其中=1n∑ni=1Xi。 2． 极大似然估计法 极大似然估计法(Maximum Likelihood Estimation)只能在已知总体分布的前提下进行，例如，假定一个盒子里装有许多大小相同的黑球和白球，并且假定它们的数目之比为3∶1，但不知是白球多还是黑球多，现在有放回地从盒中抽了3个球，试根据所抽3个球中黑球的数目确定是白球多还是黑球多． 解： 设所抽3个球中黑球数为X，摸到黑球的概率为p，则X服从二项分布 P{X=k}=Ck3pk(1-p)3-kk=0，1，2，3。 问题是p=1/4还是p=3/4呢？现根据样本中黑球数，对未知参数p进行估计。抽样后，共有4种可能结果，其概率如表31所示。 表31概论分布律 X 0 1 2 3 p=1/4时，P{X=k} 27/64 27/64 9/64 1/64 p=3/4时，P{X=k} 1/64 9/64 27/64 27/64 假如某次抽样中，只出现一个黑球，即X=1，p=1/4时，P{X=1}=27/64； p=3/4时，P{X=1}=9/64，这时我们就会选择p=1/4，即黑球数比白球数为1∶3。因为在一次试验中，事件“1个黑球”发生了。我们认为它应有较大的概率27/64(27/64＞9/64)，而27/64对应着参数p=1/4，同样可以考虑X=0，2，3的情形，最后可得 p=14当x=0，1时 34当x=2，3时 1) 似然函数 在极大似然估计法中，最关键的问题是如何求得似然函数，有了似然函数，问题就简单了，下面分两种情形来介绍似然函数。 (1) 离散型总体 设总体X为离散型，P{X=x}=p(x，θ)，其中θ为待估计的未知参数，假定x1，x2，…，xn为样本X1，X2，…，Xn的一组观测值。 P{X1=x1，X2=x2，…，Xn=xn}=P{X1=x1}P{X2=x2}…P{Xn=xn} =p(x1，θ)p(x2，θ)…p(xn，θ) =∏ni=1p(xi,θ) 将∏ni=1p(xi,θ)看作参数θ的函数，记为L(θ)，即 L(θ)=∏ni=1p(xi,θ)(312) (2) 连续型总体 设总体X为连续型，已知其分布密度函数为f(x，θ)，θ为待估计的未知参数，则样本(X1，X2，…，Xn)的联合密度为 f(x1，θ)f(x2，θ)…f(xn，θ)=∏ni=1p(xi,θ) 将它也看作关于参数θ的函数，记为L(θ)，即 L(θ)=∏ni=1p(xi,θ)(313) 由此可见，不管是离散型总体，还是连续型总体，只要知道它的概率分布或密度函数，我们总可以得到一个关于参数θ的函数L(θ)，称L(θ)为似然函数。 2) 极大似然估计 极大似然估计法的主要思想是： 如果随机抽样得到的样本观测值为x1，x2，…，xn，则我们应当这样来选取未知参数θ的值，使得出现该样本值的可能性最大，即使似然函数L(θ)取最大值，从而求参数θ的极大似然估计的问题就转化为求似然函数L(θ)的极值点的问题，一般来说，这个问题可以通过求解下面的方程来解决： dL(θ)dθ=0(314) 然而，L(θ)是n个函数的连乘积，求导数比较复杂，由于ln L(θ)是L(θ)的单调增函数，所以L(θ)与ln L(θ)在θ的同一点处取得极大值。于是求解式(33)可转化为求解： dlnL(θ)dθ=0(315) 称ln L(θ)为对数似然函数，方程(315)为对数似然方程，求解此方程就可得到参数θ的估计值。 如果总体X的分布中含有k个未知参数： θ1，θ2，…，θk，则极大似然估计法也适用。此时，所得的似然函数是关于θ1，θ2，…，θk的多元函数L(θ1，θ2，…，θk)，解下列方程组，就可得到θ1，θ2，…，θk的估计值 lnL(θ1,θ2,…,θk)θ1=0 lnL(θ1,θ2,…,θk)θ2=0  lnL(θ1,θ2,…,θk)θk=0(316) 3．3．2评价估计量的标准 设总体X服从［0，θ］上的均匀分布，θ^矩=2，θ^Lmax1≤i≤nXi都是θ的估计，这两个估计哪一个好呢？下面我们首先讨论衡量估计量好坏的标准问题。 1． 无偏性 若估计量(X1，X2，…，Xn)的数学期望等于未知参数θ，即 E(θ^)=θ(317) 则称θ^为θ的无偏估计量(Nondeviation Estimator)。 估计量θ^的值不一定就是θ的真值，因为它是一个随机变量，若θ^是θ的无偏估计，则尽管θ^的值随样本值的不同而变化，但平均来说它会等于θ的真值。 2． 有效性 对于未知参数θ，如果有两个无偏估计量θ^1与θ^2，即E(θ^1)=E(θ^2)=θ，那么在θ^1和θ^2中谁更好呢？此时我们自然希望θ的平均偏差E(θ^-θ)2越小越好，即一个好的估计量应该有尽可能小的方差，这就是有效性。 设θ^1和θ^2都是未知参数θ的无偏估计，若对任意的参数θ，有 D(θ^1)≤D(θ^2)(318) 则称θ^1比θ^2有效。 如果θ^1比θ^2有效，则虽然θ^1还不是θ的真值，但θ^1在θ附近取值的密集程度较θ^2高，即用θ^1估计θ精度要高些。 例如，对正态总体N(μ，σ2)，=1n∑ni=1Xi，Xi和都是E(X)=μ的无偏估计量，但D(X)=σ2n≤D(Xi)=σ2，故较个别观测值Xi有效。实际中也是如此，例如要估计某个班学生的平均成绩，可使用两种方法，一种方法是在该班任意抽取一个同学，以该同学的成绩作为全班的平均成绩； 另一种方法是在该班抽取n位同学，以这n个同学的平均成绩作为全班的平均成绩，显然第二种方法比第一种方法好。 3． 一致性 无偏性、有效性都是在样本容量n一定的条件下进行讨论的，然而(X1，X2，…，Xn)不仅与样本值有关，而且与样本容量n有关，记为n，很自然，我们希望n越大时，n对θ的估计应该越精确。如果n依概率收敛于θ，即ε＞0，有 limn→∞Pθ^n－θ<ε=1(319) 则称θ^n是θ的一致估计量(Uniform Estimator)。 3．3．3区间估计 1． 区间估计的概念 3.3.1节介绍了参数的点估计，假设总体X~N(μ，σ2)，对于样本(X1，X2，…，Xn)，μ^=是参数μ的矩法估计和极大似然估计，并且满足无偏性和一致性。但实际上=μ的可能性有多大呢？由于是一连续型随机变量，P{X=μ}=0，即μ^=μ的可能性为0，因此，我们希望给出μ的一个大致范围，使得μ有较高的概率在这个范围内，这就是区间估计问题。 设θ^1(X1，X2，…，Xn)及θ^2 (X1，X2，…，Xn)是两个统计量，如果对于给定的概率1-α(0<α<1)，有 P{θ^1＜θ＜θ^2}=1-α(320) 则称随机区间(θ^1，θ^2)为参数θ的置信区间(Confidence Interval)，θ^1称为置信下限，θ^2称为置信上限，1-α称为置信概率或置信度(Confidence Level)。 定义中的随机区间(θ^1，θ^2)的大小依赖于随机抽取的样本观测值，它可能包含θ，也可能不包含θ，式(39)的意义是指(θ^1，θ^2)以1-α的概率包含θ。例如，若取α=0．05，那么置信概率为1-α=0．95，这时，置信区间(θ^1，θ^2)的意义是指： 在100次重复抽样所得到的100个置信区间中，大约有95个区间包含参数真值θ，有5个区间不包含真值θ，亦即随机区间(θ^1，θ^2)包含参数θ真值的频率近似为0．95。 2． 正态总体参数的区间估计 由于在大多数情况下，所遇到的总体是服从正态分布的(有的是近似正态分布)，下面重点讨论正态总体参数的区间估计问题。在下面的讨论中，总假定X~N(μ，σ2)，X1，X2，…，Xn为其样本。 分两种情况进行讨论。如果σ2已知，则μ的置信区间为－zα2σn，+zα2σn， 置信概率为1-α。 如果σ2未知，不能使用式(37)作为置信区间，因为式(37)中区间的端点与σ有关，考虑到S2=1n－1∑ni=1(Xi－)2是σ2的无偏估计，将－μσ/n中的σ换成S得 T=－μS/n ~t(n-1) 对于给定的α，由t分布分位数表可得上分位点tσ/2(n-1)，使得 P－μS/n