statistics
— 数学统计函数¶
3.4 版新增。
源代码: Lib/statistics.py
此模块提供用于计算数值(Real
值)数据的数学统计量的函数。
该模块并非旨在与第三方库(如 NumPy、SciPy)或面向专业统计学家的专有全功能统计软件包(如 Minitab、SAS 和 Matlab)竞争。它的目标是图形计算器和科学计算器的水平。
除非另有说明,否则这些函数支持 int
、float
、Decimal
和 Fraction
。目前不支持对其他类型(无论是否在数字塔中)的行为。混合类型的集合也是未定义的,并且取决于实现。如果您的输入数据由混合类型组成,则可以使用 map()
来确保一致的结果,例如:map(float, input_data)
。
某些数据集使用 NaN
(非数字)值来表示缺失数据。由于 NaN 具有不寻常的比较语义,因此它们会在对数据进行排序或计算出现次数的统计函数中导致意外或未定义的行为。受影响的函数有 median()
、median_low()
、median_high()
、median_grouped()
、mode()
、multimode()
和 quantiles()
。在调用这些函数之前,应删除 NaN
值。
>>> from statistics import median
>>> from math import isnan
>>> from itertools import filterfalse
>>> data = [20.7, float('NaN'),19.2, 18.3, float('NaN'), 14.4]
>>> sorted(data) # This has surprising behavior
[20.7, nan, 14.4, 18.3, 19.2, nan]
>>> median(data) # This result is unexpected
16.35
>>> sum(map(isnan, data)) # Number of missing values
2
>>> clean = list(filterfalse(isnan, data)) # Strip NaN values
>>> clean
[20.7, 19.2, 18.3, 14.4]
>>> sorted(clean) # Sorting now works as expected
[14.4, 18.3, 19.2, 20.7]
>>> median(clean) # This result is now well defined
18.75
平均值和集中趋势度量¶
这些函数从总体或样本中计算平均值或典型值。
数据的算术平均值(“平均数”)。 |
|
快速浮点算术平均值,支持可选的加权。 |
|
数据的几何平均值。 |
|
数据的调和平均值。 |
|
数据的中位数(中间值)。 |
|
数据的低位中位数。 |
|
数据的高位中位数。 |
|
分组数据的中间值(第 50 个百分位数)。 |
|
离散或名义数据的单一众数(最常见的值)。 |
|
离散或名义数据的众数列表(最常见的值)。 |
|
将数据分成概率相等的区间。 |
离散度度量¶
这些函数计算总体或样本偏离典型值或平均值的程度。
数据的总体标准差。 |
|
数据的总体方差。 |
|
数据的样本标准差。 |
|
数据的样本方差。 |
两个输入之间关系的统计量¶
这些函数计算关于两个输入之间关系的统计量。
两个变量的样本协方差。 |
|
皮尔逊和斯皮尔曼相关系数。 |
|
简单线性回归的斜率和截距。 |
函数细节¶
注意:这些函数不要求提供给它们的数据已排序。但是,为了便于阅读,大多数示例都显示了已排序的序列。
- statistics.mean(data)¶
返回 data 的样本算术平均值,它可以是序列或可迭代对象。
算术平均值是数据的总和除以数据点的数量。它通常被称为“平均数”,尽管它只是许多不同数学平均值中的一种。它是数据集中趋势的度量。
如果 data 为空,则会引发
StatisticsError
。一些使用示例
>>> mean([1, 2, 3, 4, 4]) 2.8 >>> mean([-1.0, 2.5, 3.25, 5.75]) 2.625 >>> from fractions import Fraction as F >>> mean([F(3, 7), F(1, 21), F(5, 3), F(1, 3)]) Fraction(13, 21) >>> from decimal import Decimal as D >>> mean([D("0.5"), D("0.75"), D("0.625"), D("0.375")]) Decimal('0.5625')
- statistics.fmean(data, weights=None)¶
将 data 转换为浮点数并计算算术平均值。
这比
mean()
函数运行速度更快,并且它始终返回float
。data 可以是序列或可迭代对象。如果输入数据集为空,则引发StatisticsError
。>>> fmean([3.5, 4.0, 5.25]) 4.25
支持可选的加权。例如,教授通过将测验加权 20%、家庭作业加权 20%、期中考试加权 30% 和期末考试加权 30% 来为课程评分。
>>> grades = [85, 92, 83, 91] >>> weights = [0.20, 0.20, 0.30, 0.30] >>> fmean(grades, weights) 87.6
如果提供了 weights,则其长度必须与 data 相同,否则将引发
ValueError
。3.8 版新增。
在 3.11 版更改: 添加了对 weights 的支持。
- statistics.geometric_mean(data)¶
将 data 转换为浮点数并计算几何平均值。
几何平均值使用值的乘积(而不是使用其总和的算术平均值)来指示 data 的集中趋势或典型值。
如果输入数据集为空、包含零或包含负值,则引发
StatisticsError
。data 可以是序列或可迭代对象。没有特别努力去实现精确的结果。(但是,这在将来可能会改变。)
>>> round(geometric_mean([54, 24, 36]), 1) 36.0
3.8 版新增。
- statistics.harmonic_mean(data, weights=None)¶
返回 data 的调和平均数,data 是一个由实值数组成的序列或可迭代对象。如果省略 weights 或为
None
,则假设权重相等。调和平均数是数据倒数的算术
mean()
的倒数。例如,三个值 a、b 和 c 的调和平均数等价于3/(1/a + 1/b + 1/c)
。如果其中一个值为零,则结果为零。调和平均数是一种平均数,是数据集中趋势的度量。它通常适用于对比率或速率(例如速度)进行平均。
假设一辆汽车以 40 公里/小时的速度行驶 10 公里,然后以 60 公里/小时的速度行驶另外 10 公里。平均速度是多少?
>>> harmonic_mean([40, 60]) 48.0
假设一辆汽车以 40 公里/小时的速度行驶 5 公里,当交通顺畅时,在剩余的 30 公里路程中加速到 60 公里/小时。平均速度是多少?
>>> harmonic_mean([40, 60], weights=[5, 30]) 56.0
如果 data 为空、任何元素小于零或加权和不为正,则会引发
StatisticsError
。当前算法在遇到输入中的零时会提前退出。这意味着不会测试后续输入的有效性。(此行为将来可能会更改。)
在 3.6 版中添加。
在 3.10 版中更改: 添加了对 weights 的支持。
- statistics.median(data)¶
使用常见的“中间两个数的平均值”方法返回数值数据的中位数(中间值)。如果 data 为空,则会引发
StatisticsError
。data 可以是序列或可迭代对象。中位数是集中趋势的可靠度量,受异常值的影响较小。当数据点的数量为奇数时,返回中间数据点
>>> median([1, 3, 5]) 3
当数据点的数量为偶数时,通过取中间两个值的平均值来插值中位数
>>> median([1, 3, 5, 7]) 4.0
这适用于数据是离散的,并且您不介意中位数可能不是实际数据点的情况。
如果数据是有序的(支持顺序操作)但不是数值的(不支持加法),请考虑改用
median_low()
或median_high()
。
- statistics.median_low(data)¶
返回数值数据的低位中位数。如果 data 为空,则会引发
StatisticsError
。data 可以是序列或可迭代对象。低位中位数始终是数据集的成员。当数据点的数量为奇数时,返回中间值。当它为偶数时,返回两个中间值中较小的一个。
>>> median_low([1, 3, 5]) 3 >>> median_low([1, 3, 5, 7]) 3
当您的数据是离散的,并且您希望中位数是实际数据点而不是插值时,请使用低位中位数。
- statistics.median_high(data)¶
返回数据的高位中位数。如果 data 为空,则会引发
StatisticsError
。data 可以是序列或可迭代对象。高位中位数始终是数据集的成员。当数据点的数量为奇数时,返回中间值。当它为偶数时,返回两个中间值中较大的一个。
>>> median_high([1, 3, 5]) 3 >>> median_high([1, 3, 5, 7]) 5
当您的数据是离散的,并且您希望中位数是实际数据点而不是插值时,请使用高位中位数。
- statistics.median_grouped(data, interval=1.0)¶
估计已围绕连续的、固定宽度的间隔的中点进行 分组或分箱 的数值数据的中位数。
data 可以是任何数值数据的可迭代对象,每个值都恰好是分箱的中点。必须至少存在一个值。
interval 是每个分箱的宽度。
例如,人口统计信息可能已汇总到连续的十年年龄组中,每个组由间隔的 5 年中点表示
>>> from collections import Counter >>> demographics = Counter({ ... 25: 172, # 20 to 30 years old ... 35: 484, # 30 to 40 years old ... 45: 387, # 40 to 50 years old ... 55: 22, # 50 to 60 years old ... 65: 6, # 60 to 70 years old ... }) ...
第 50 个百分位数(中位数)是 1071 名成员队列中的第 536 个人。该人属于 30 至 40 岁的年龄组。
常规的
median()
函数会假设三十多岁年龄组中的每个人都正好 35 岁。更合理的假设是,该年龄组的 484 名成员均匀分布在 30 到 40 岁之间。为此,我们使用median_grouped()
>>> data = list(demographics.elements()) >>> median(data) 35 >>> round(median_grouped(data, interval=10), 1) 37.5
调用者负责确保数据点之间以 interval 的精确倍数分隔。这对于获得正确的结果至关重要。该函数不检查此先决条件。
输入可以是在插值步骤中可以强制转换为浮点数的任何数值类型。
- statistics.mode(data)¶
从离散或名义 data 中返回单个最常见的数据点。众数(如果存在)是最典型的值,用作集中趋势的度量。
如果有多个具有相同频率的众数,则返回在 data 中遇到的第一个众数。如果需要最小或最大的那个,请使用
min(multimode(data))
或max(multimode(data))
。如果输入 data 为空,则会引发StatisticsError
。mode
假设数据是离散的,并返回单个值。这是学校通常教授的众数的标准处理方式>>> mode([1, 1, 2, 3, 3, 3, 3, 4]) 3
众数的独特之处在于它是此包中唯一也适用于名义(非数值)数据的统计量
>>> mode(["red", "blue", "blue", "red", "green", "red", "red"]) 'red'
在 3.8 版中更改: 现在通过返回遇到的第一个众数来处理多峰数据集。以前,当找到多个众数时,它会引发
StatisticsError
。
- statistics.multimode(data)¶
按它们在 data 中首次出现的顺序返回最常出现的值的列表。如果有多个众数,则返回多个结果;如果 data 为空,则返回空列表
>>> multimode('aabbbbccddddeeffffgg') ['b', 'd', 'f'] >>> multimode('') []
3.8 版新增。
- statistics.pstdev(data, mu=None)¶
返回总体标准差(总体方差的平方根)。有关参数和其他详细信息,请参阅
pvariance()
。>>> pstdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75]) 0.986893273527251
- statistics.pvariance(data, mu=None)¶
返回 *data* 的总体方差,*data* 是一个非空的实值数字序列或可迭代对象。方差,或关于均值的二阶矩,是衡量数据可变性(散布或离散程度)的指标。大方差表示数据分散;小方差表示数据紧密地聚集在均值周围。
如果给出了可选的第二个参数 *mu*,则它应该是 *data* 的 *总体* 均值。它也可以用来计算围绕非均值点的二阶矩。如果缺少或为
None
(默认值),则自动计算算术平均值。使用此函数可以根据整个总体计算方差。要根据样本估计方差,
variance()
函数通常是更好的选择。如果 *data* 为空,则引发
StatisticsError
。示例
>>> data = [0.0, 0.25, 0.25, 1.25, 1.5, 1.75, 2.75, 3.25] >>> pvariance(data) 1.25
如果您已经计算了数据的均值,则可以将其作为可选的第二个参数 *mu* 传递,以避免重新计算
>>> mu = mean(data) >>> pvariance(data, mu) 1.25
支持小数和分数
>>> from decimal import Decimal as D >>> pvariance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")]) Decimal('24.815') >>> from fractions import Fraction as F >>> pvariance([F(1, 4), F(5, 4), F(1, 2)]) Fraction(13, 72)
注意
当使用整个总体调用时,这将给出总体方差 σ²。当改为对样本调用时,这是有偏样本方差 s²,也称为具有 N 个自由度的方差。
如果您以某种方式知道真实的总体均值 μ,则可以使用此函数来计算样本的方差,将已知的总体均值作为第二个参数给出。如果数据点是总体的随机样本,则结果将是总体方差的无偏估计。
- statistics.stdev(data, xbar=None)¶
返回样本标准差(样本方差的平方根)。有关参数和其他详细信息,请参阅
variance()
。>>> stdev([1.5, 2.5, 2.5, 2.75, 3.25, 4.75]) 1.0810874155219827
- statistics.variance(data, xbar=None)¶
返回 *data* 的样本方差,*data* 是一个至少包含两个实值数字的可迭代对象。方差,或关于均值的二阶矩,是衡量数据可变性(散布或离散程度)的指标。大方差表示数据分散;小方差表示数据紧密地聚集在均值周围。
如果给出了可选的第二个参数 *xbar*,则它应该是 *data* 的 *样本* 均值。如果缺少或为
None
(默认值),则自动计算均值。当您的数据是总体中的样本时,请使用此函数。要根据整个总体计算方差,请参阅
pvariance()
。如果 *data* 的值少于两个,则引发
StatisticsError
。示例
>>> data = [2.75, 1.75, 1.25, 0.25, 0.5, 1.25, 3.5] >>> variance(data) 1.3720238095238095
如果您已经计算了数据的样本均值,则可以将其作为可选的第二个参数 *xbar* 传递,以避免重新计算
>>> m = mean(data) >>> variance(data, m) 1.3720238095238095
此函数不会尝试验证您是否已将实际均值作为 *xbar* 传递。对 *xbar* 使用任意值可能会导致无效或不可能的结果。
支持小数和分数
>>> from decimal import Decimal as D >>> variance([D("27.5"), D("30.25"), D("30.25"), D("34.5"), D("41.75")]) Decimal('31.01875') >>> from fractions import Fraction as F >>> variance([F(1, 6), F(1, 2), F(5, 3)]) Fraction(67, 108)
注意
这是具有贝塞尔校正的样本方差 s²,也称为具有 N-1 个自由度的方差。如果数据点具有代表性(例如,独立且同分布),则结果应该是真实总体方差的无偏估计。
如果您以某种方式知道实际总体均值 μ,则应将其作为 *mu* 参数传递给
pvariance()
函数,以获取样本的方差。
- statistics.quantiles(data, *, n=4, method='exclusive')¶
将 *data* 分成 *n* 个概率相等的连续区间。返回一个包含
n - 1
个分割点的列表,这些分割点用于分隔区间。将 *n* 设置为 4 表示四分位数(默认值)。将 *n* 设置为 10 表示十分位数。将 *n* 设置为 100 表示百分位数,这将给出 99 个分割点,将 *data* 分成 100 个大小相等的组。如果 *n* 不小于 1,则引发
StatisticsError
。*data* 可以是任何包含样本数据的可迭代对象。为了获得有意义的结果,*data* 中的数据点数量应大于 *n*。如果数据点少于两个,则引发
StatisticsError
。分割点由两个最近的数据点线性插值得到。例如,如果一个分割点位于两个样本值(
100
和112
)之间距离的三分之一处,则该分割点的值为104
。计算分位数的 *method* 可以有所不同,具体取决于 *data* 是否包含总体中可能存在的最小值和最大值。
默认的 *method* 是 “exclusive”,用于从总体中采样的数据,这些数据的极值可能比样本中的极值更极端。低于 *m* 个排序数据点中第 *i* 个数据点的总体部分计算为
i / (m + 1)
。给定九个样本值,该方法对它们进行排序并分配以下百分位数:10%、20%、30%、40%、50%、60%、70%、80%、90%。将method设置为“inclusive”用于描述总体数据或已知包含总体中最值样本的情况。data中的最小值被视为第 0 个百分位数,最大值被视为第 100 个百分位数。低于m个排序数据点中第i个的总体部分计算为
(i - 1) / (m - 1)
。给定 11 个样本值,该方法将对它们进行排序并分配以下百分位数:0%、10%、20%、30%、40%、50%、60%、70%、80%、90%、100%。# Decile cut points for empirically sampled data >>> data = [105, 129, 87, 86, 111, 111, 89, 81, 108, 92, 110, ... 100, 75, 105, 103, 109, 76, 119, 99, 91, 103, 129, ... 106, 101, 84, 111, 74, 87, 86, 103, 103, 106, 86, ... 111, 75, 87, 102, 121, 111, 88, 89, 101, 106, 95, ... 103, 107, 101, 81, 109, 104] >>> [round(q, 1) for q in quantiles(data, n=10)] [81.0, 86.2, 89.0, 99.4, 102.5, 103.6, 106.0, 109.8, 111.0]
3.8 版新增。
- statistics.covariance(x, y, /)¶
返回两个输入x和y的样本协方差。协方差是两个输入联合变异性的度量。
两个输入的长度必须相同(不少于两个),否则将引发
StatisticsError
。示例
>>> x = [1, 2, 3, 4, 5, 6, 7, 8, 9] >>> y = [1, 2, 3, 1, 2, 3, 1, 2, 3] >>> covariance(x, y) 0.75 >>> z = [9, 8, 7, 6, 5, 4, 3, 2, 1] >>> covariance(x, z) -7.5 >>> covariance(z, x) -7.5
版本 3.10 中的新功能。
- statistics.correlation(x, y, /, *, method='linear')¶
返回两个输入的皮尔逊相关系数。皮尔逊相关系数r的取值范围为 -1 到 +1。它衡量线性关系的强度和方向。
如果method为“ranked”,则计算两个输入的斯皮尔曼等级相关系数。数据将被等级替换。并列排名将被平均,以便相等的值获得相同的排名。所得系数衡量单调关系的强度。
斯皮尔曼相关系数适用于序数数据或不满足皮尔逊相关系数线性比例要求的连续数据。
两个输入的长度必须相同(不少于两个),并且不能为常数,否则将引发
StatisticsError
。>>> # Mercury, Venus, Earth, Mars, Jupiter, Saturn, Uranus, and Neptune >>> orbital_period = [88, 225, 365, 687, 4331, 10_756, 30_687, 60_190] # days >>> dist_from_sun = [58, 108, 150, 228, 778, 1_400, 2_900, 4_500] # million km >>> # Show that a perfect monotonic relationship exists >>> correlation(orbital_period, dist_from_sun, method='ranked') 1.0 >>> # Observe that a linear relationship is imperfect >>> round(correlation(orbital_period, dist_from_sun), 4) 0.9882 >>> # Demonstrate Kepler's third law: There is a linear correlation >>> # between the square of the orbital period and the cube of the >>> # distance from the sun. >>> period_squared = [p * p for p in orbital_period] >>> dist_cubed = [d * d * d for d in dist_from_sun] >>> round(correlation(period_squared, dist_cubed), 4) 1.0
版本 3.10 中的新功能。
版本 3.12 中的变化: 添加了对斯皮尔曼等级相关系数的支持。
- statistics.linear_regression(x, y, /, *, proportional=False)¶
返回使用普通最小二乘法估计的简单线性回归参数的斜率和截距。简单线性回归用以下线性函数描述自变量x和因变量y之间的关系
y = slope * x + intercept + noise
其中
slope
和intercept
是估计的回归参数,noise
表示未被线性回归解释的数据的变异性(它等于因变量的预测值与实际值之间的差)。两个输入的长度必须相同(不少于两个),并且自变量x不能为常数;否则将引发
StatisticsError
。例如,我们可以使用巨蟒电影的上映日期来预测到 2019 年巨蟒电影的累计产量,假设它们保持这种速度。
>>> year = [1971, 1975, 1979, 1982, 1983] >>> films_total = [1, 2, 3, 4, 5] >>> slope, intercept = linear_regression(year, films_total) >>> round(slope * 2019 + intercept) 16
如果proportional为 true,则假设自变量x和因变量y成正比。数据拟合到一条穿过原点的线。由于intercept始终为 0.0,因此基础线性函数简化为
y = slope * x + noise
继续
correlation()
中的示例,我们来看看基于主要行星的模型在预测矮行星的轨道距离方面的效果如何>>> model = linear_regression(period_squared, dist_cubed, proportional=True) >>> slope = model.slope >>> # Dwarf planets: Pluto, Eris, Makemake, Haumea, Ceres >>> orbital_periods = [90_560, 204_199, 111_845, 103_410, 1_680] # days >>> predicted_dist = [math.cbrt(slope * (p * p)) for p in orbital_periods] >>> list(map(round, predicted_dist)) [5912, 10166, 6806, 6459, 414] >>> [5_906, 10_152, 6_796, 6_450, 414] # actual distance in million km [5906, 10152, 6796, 6450, 414]
版本 3.10 中的新功能。
版本 3.11 中的变化: 添加了对proportional的支持。
异常¶
定义了一个异常
- exception statistics.StatisticsError¶
ValueError
的子类,用于与统计相关的异常。
NormalDist
对象¶
NormalDist
是一种用于创建和操作随机变量的正态分布的工具。它是一个将数据测量的均值和标准差视为单个实体的类。
正态分布源于中心极限定理,在统计学中有着广泛的应用。
- class statistics.NormalDist(mu=0.0, sigma=1.0)¶
返回一个新的NormalDist对象,其中mu表示算术平均值,sigma表示标准差。
如果sigma为负数,则引发
StatisticsError
。- classmethod from_samples(data)¶
使用
fmean()
和stdev()
从 *data* 估计 *mu* 和 *sigma* 参数,创建一个正态分布实例。*data* 可以是任何可迭代对象,并且应该包含可以转换为
float
类型的数值。如果 *data* 不包含至少两个元素,则会引发StatisticsError
,因为估计中心值至少需要一个点,估计离散度至少需要两个点。
- samples(n, *, seed=None)¶
为给定的均值和标准差生成 *n* 个随机样本。返回一个包含
float
值的list
。如果给定 *seed*,则创建一个新的底层随机数生成器实例。这对于创建可重复的结果非常有用,即使在多线程环境中也是如此。
- pdf(x)¶
使用概率密度函数 (pdf),计算随机变量 *X* 接近给定值 *x* 的相对可能性。在数学上,它是当 *dx* 趋近于零时,比率
P(x <= X < x+dx) / dx
的极限。相对可能性计算为样本出现在一个狭窄范围内的概率除以范围的宽度(因此称为“密度”)。由于可能性是相对于其他点的,因此其值可能大于
1.0
。
- cdf(x)¶
使用累积分布函数 (cdf),计算随机变量 *X* 小于或等于 *x* 的概率。在数学上,它写成
P(X <= x)
。
- inv_cdf(p)¶
计算逆累积分布函数,也称为分位数函数或百分点函数。在数学上,它写成
x : P(X <= x) = p
。找到随机变量 *X* 的值 *x*,使得变量小于或等于该值的概率等于给定的概率 *p*。
- overlap(other)¶
测量两个正态概率分布之间的一致性。返回一个介于 0.0 和 1.0 之间的值,表示两个概率密度函数的重叠区域。
- quantiles(n=4)¶
将正态分布分成 *n* 个概率相等的连续区间。返回一个包含 (n - 1) 个分割点的列表,这些分割点将区间分开。
将 *n* 设置为 4 表示四分位数(默认值)。将 *n* 设置为 10 表示十分位数。将 *n* 设置为 100 表示百分位数,这将给出 99 个分割点,将正态分布分成 100 个大小相等的组。
NormalDist
的实例支持与常数的加法、减法、乘法和除法运算。这些运算用于平移和缩放。例如>>> temperature_february = NormalDist(5, 2.5) # Celsius >>> temperature_february * (9/5) + 32 # Fahrenheit NormalDist(mu=41.0, sigma=4.5)
不支持用常数除以
NormalDist
的实例,因为结果将不是正态分布的。由于正态分布是由自变量的累加效应产生的,因此可以对表示为
NormalDist
实例的两个独立正态分布随机变量进行加减运算。例如>>> birth_weights = NormalDist.from_samples([2.5, 3.1, 2.1, 2.4, 2.7, 3.5]) >>> drug_effects = NormalDist(0.4, 0.15) >>> combined = birth_weights + drug_effects >>> round(combined.mean, 1) 3.1 >>> round(combined.stdev, 1) 0.5
3.8 版新增。
示例和方法¶
经典概率问题¶
NormalDist
可以轻松解决经典概率问题。
例如,根据SAT 考试的历史数据显示,分数呈正态分布,均值为 1060,标准差为 195,确定分数在 1100 到 1200 之间的学生所占的百分比,四舍五入到最接近的整数
>>> sat = NormalDist(1060, 195)
>>> fraction = sat.cdf(1200 + 0.5) - sat.cdf(1100 - 0.5)
>>> round(fraction * 100.0, 1)
18.4
>>> list(map(round, sat.quantiles()))
[928, 1060, 1192]
>>> list(map(round, sat.quantiles(n=10)))
[810, 896, 958, 1011, 1060, 1109, 1162, 1224, 1310]
模拟的蒙特卡洛输入¶
为了估计难以解析求解的模型的分布,NormalDist
可以为 蒙特卡洛模拟 生成输入样本。
>>> def model(x, y, z):
... return (3*x + 7*x*y - 5*y) / (11 * z)
...
>>> n = 100_000
>>> X = NormalDist(10, 2.5).samples(n, seed=3652260728)
>>> Y = NormalDist(15, 1.75).samples(n, seed=4582495471)
>>> Z = NormalDist(50, 1.25).samples(n, seed=6582483453)
>>> quantiles(map(model, X, Y, Z))
[1.4591308524824727, 1.8035946855390597, 2.175091447274739]
近似二项分布¶
当样本量较大且成功试验的概率接近 50% 时,可以使用正态分布来近似 二项分布。
例如,一个开源会议有 750 名与会者和两个可容纳 500 人的房间。有一个关于 Python 的演讲和一个关于 Ruby 的演讲。在之前的会议中,65% 的与会者更喜欢听 Python 演讲。假设人群偏好没有改变,那么 Python 房间保持在其容量限制内的概率是多少?
>>> n = 750 # Sample size
>>> p = 0.65 # Preference for Python
>>> q = 1.0 - p # Preference for Ruby
>>> k = 500 # Room capacity
>>> # Approximation using the cumulative normal distribution
>>> from math import sqrt
>>> round(NormalDist(mu=n*p, sigma=sqrt(n*p*q)).cdf(k + 0.5), 4)
0.8402
>>> # Exact solution using the cumulative binomial distribution
>>> from math import comb, fsum
>>> round(fsum(comb(n, r) * p**r * q**(n-r) for r in range(k+1)), 4)
0.8402
>>> # Approximation using a simulation
>>> from random import seed, binomialvariate
>>> seed(8675309)
>>> mean(binomialvariate(n, p) <= k for i in range(10_000))
0.8406
朴素贝叶斯分类器¶
正态分布在机器学习问题中很常见。
维基百科有一个关于 朴素贝叶斯分类器 的很好的例子。挑战在于从包括身高、体重和脚尺寸在内的正态分布特征的测量值来预测一个人的性别。
我们有一个包含八个人测量的训练数据集。假设测量值呈正态分布,因此我们使用 NormalDist
来汇总数据。
>>> height_male = NormalDist.from_samples([6, 5.92, 5.58, 5.92])
>>> height_female = NormalDist.from_samples([5, 5.5, 5.42, 5.75])
>>> weight_male = NormalDist.from_samples([180, 190, 170, 165])
>>> weight_female = NormalDist.from_samples([100, 150, 130, 150])
>>> foot_size_male = NormalDist.from_samples([12, 11, 12, 10])
>>> foot_size_female = NormalDist.from_samples([6, 8, 7, 9])
接下来,我们遇到一个新的人,他的特征测量值是已知的,但性别未知。
>>> ht = 6.0 # height
>>> wt = 130 # weight
>>> fs = 8 # foot size
从男性或女性的 50% 先验概率 开始,我们将后验概率计算为先验概率乘以给定性别的特征测量的似然积。
>>> prior_male = 0.5
>>> prior_female = 0.5
>>> posterior_male = (prior_male * height_male.pdf(ht) *
... weight_male.pdf(wt) * foot_size_male.pdf(fs))
>>> posterior_female = (prior_female * height_female.pdf(ht) *
... weight_female.pdf(wt) * foot_size_female.pdf(fs))
最终预测指向最大的后验概率。这被称为 最大后验概率 或 MAP。
>>> 'male' if posterior_male > posterior_female else 'female'
'female'
核密度估计¶
可以从固定数量的离散样本估计连续概率分布。
基本思想是使用 核函数(如正态分布、三角分布或均匀分布)对数据进行平滑处理。平滑程度由缩放参数 h
控制,该参数称为*带宽*。
from random import choice, random
def kde_normal(data, h):
"Create a continuous probability distribution from discrete samples."
# Smooth the data with a normal distribution kernel scaled by h.
K_h = NormalDist(0.0, h)
def pdf(x):
'Probability density function. P(x <= X < x+dx) / dx'
return sum(K_h.pdf(x - x_i) for x_i in data) / len(data)
def cdf(x):
'Cumulative distribution function. P(X <= x)'
return sum(K_h.cdf(x - x_i) for x_i in data) / len(data)
def rand():
'Random selection from the probability distribution.'
return choice(data) + K_h.inv_cdf(random())
return pdf, cdf, rand
维基百科有一个例子,我们可以使用 kde_normal()
配方来生成和绘制从一个小样本估计的概率密度函数。
>>> sample = [-2.1, -1.3, -0.4, 1.9, 5.1, 6.2]
>>> pdf, cdf, rand = kde_normal(sample, h=1.5)
>>> xarr = [i/100 for i in range(-750, 1100)]
>>> yarr = [pdf(x) for x in xarr]
xarr
和 yarr
中的点可用于绘制 PDF 图。
重采样 数据以生成 100 个新选择。
>>> new_selections = [rand() for i in range(100)]
确定新选择低于 2.0
的概率。
>>> round(cdf(2.0), 4)
0.5794
添加一个新的样本数据点并找到 2.0
处的新 CDF。
>>> sample.append(4.9)
>>> round(cdf(2.0), 4)
0.5005