tnkstat3e-merge-1

统计思维（程序员的概率统计）第三版（二）

原文：allendowney.github.io/ThinkStats/index.html

译者：飞龙

协议：CC BY-NC-SA 4.0

模型分布

原文：allendowney.github.io/ThinkStats/chap05.html

我们迄今为止使用的分布被称为经验分布，因为它们基于经验观察——换句话说，数据。我们在现实世界中看到的许多数据集都可以被一个理论分布紧密近似，这通常基于一个简单的数学函数。本章介绍了这些理论分布及其可以用来建模的数据集。

作为例子，我们将看到：

在飞碟射击比赛中，命中和未命中的次数可以用二项分布很好地描述。
在像曲棍球和足球（足球）这样的游戏中，一场比赛中的进球数遵循泊松分布，进球之间的时间间隔遵循指数分布。
出生体重遵循正态分布，也称为高斯分布，而成年体重遵循对数正态分布。

如果你对这些分布——或者这些运动——不熟悉，我将解释你需要知道的内容。对于每个例子，我们将从一个基于简单模型的模拟开始，并展示模拟结果遵循理论分布。然后我们将看到实际数据与模型的一致性如何。

点击此处运行此笔记本在 Colab。

from  os.path  import basename, exists

def  download(url):
    filename = basename(url)
    if not exists(filename):
        from  urllib.request  import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("Downloaded " + local)

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")

try:
    import  empiricaldist
except ImportError:
    %pip install empiricaldist

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  plt

from  thinkstats  import decorate

二项分布

作为第一个例子，我们将考虑飞碟射击这项运动，其中参赛者使用猎枪射击被抛入空中的粘土碟。在国际比赛中，包括奥运会，有五轮，每轮 25 个目标，根据需要增加轮次以确定胜者。

作为飞碟射击比赛的模型，假设每个参与者命中每个目标的概率相同，p。当然，这个模型是一种简化——在现实中，一些参赛者的概率比其他人高，即使是单个参赛者，也可能从一次尝试到下一次尝试而变化。但即使它不现实，这个模型也能做出一些令人惊讶的准确预测，正如我们将看到的。

为了模拟这个模型，我将使用以下函数，它接受目标数n和每个目标的命中概率p，并返回一个由 1 和 0 组成的序列，以表示命中和未命中。

def  flip(n, p):
    choices = [1, 0]
    probs = [p, 1 - p]
    return np.random.choice(choices, n, p=probs)

这里有一个模拟 25 个目标的例子，其中每个目标的命中概率为 90%。

# Seed the random number generator so we get the same results every time
np.random.seed(1)

flip(25, 0.9)

array([1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0,
       1, 1, 1])

如果我们生成一个更长的序列并计算结果的Pmf，我们可以确认 1 和 0 的比例是正确的，至少是近似正确的。

from  empiricaldist  import Pmf

seq = flip(1000, 0.9)
pmf = Pmf.from_seq(seq)
pmf

	probs
0	0.101
1	0.899

现在我们可以使用flip来模拟一轮飞碟射击并返回命中的次数。

def  simulate_round(n, p):
    seq = flip(n, p)
    return seq.sum()

在一场大型比赛中，假设 200 名参赛者每人射击 5 轮，所有射击命中目标的概率相同，p=0.9。我们可以通过调用simulate_round 1000 次来模拟这样的比赛。

n = 25
p = 0.9
results_sim = [simulate_round(n, p) for i in range(1000)]

平均分接近22.5，这是n和p的乘积。

np.mean(results_sim), n * p

(np.float64(22.522), 22.5)

下面是结果分布的图示。

from  empiricaldist  import Pmf

pmf_sim = Pmf.from_seq(results_sim, name="simulation results")

pmf_sim.bar()
decorate(xlabel="Hits", ylabel="PMF")

顶峰接近平均值，分布向左偏斜。

我们可以不运行模拟，而是预测这个分布。从数学上讲，这些结果的分布遵循二项分布，其 PMF 易于计算。

from  scipy.special  import comb

def  binomial_pmf(k, n, p):
    return comb(n, k) * (p**k) * ((1 - p) ** (n - k))

SciPy 提供了comb函数，该函数计算从n个物品中每次取k个的组合数，通常读作“n 选 k”。

binomial_pmf计算在给定p的情况下，从n次尝试中得到k次命中的概率。如果我们用一系列的k值调用此函数，我们可以创建一个表示结果分布的Pmf。

ks = np.arange(16, n + 1)
ps = binomial_pmf(ks, n, p)
pmf_binom = Pmf(ps, ks, name="binomial model")

下面是它与模拟结果相比的图示。

from  thinkstats  import two_bar_plots

two_bar_plots(pmf_sim, pmf_binom)
decorate(xlabel="Hits", ylabel="PMF")

它们很相似，由于模拟结果中的随机变化，存在一些小的差异。这种一致性并不令人惊讶，因为模拟和模型基于相同的假设——尤其是每个尝试成功的概率相同的假设。对模型的一个更强测试是它与实际数据的比较。

从 2020 年夏季奥运会男子飞碟射击比赛的维基百科页面，我们可以提取一个表格，显示资格赛的结果。下载数据的说明在本书的笔记本中。

从en.wikipedia.org/wiki/Shooting_at_the_2020_Summer_Olympics_%E2%80%93_Men's_skeet于 2024 年 7 月 15 日下载。

filename = "Shooting_at_the_2020_Summer_Olympics_Mens_skeet"

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/" + filename)

tables = pd.read_html(filename)
table = tables[6]
table.head()

	Rank	Athlete	Country	1	2	3	4	5	Total[3]	Shoot-off	Notes
0	1	Éric Delaunay	France	25	25	25	24	25	124	+6	Q, OR
1	2	Tammaro Cassandro	Italy	24	25	25	25	25	124	+5	Q, OR
2	3	Eetu Kallioinen	Finland	25	25	24	25	24	123	NaN	Q
3	4	Vincent Hancock	United States	25	25	25	25	22	122	+8	Q
4	5	Abdullah Al-Rashidi	Kuwait	25	25	24	25	23	122	+7	Q

表格为每位选手一行，每轮比赛有一个列。我们将选择包含这些结果的列，并使用 NumPy 函数flatten将它们放入一个单一的数组中。

columns = ["1", "2", "3", "4", "5"]
results = table[columns].values.flatten()

有 30 名选手，我们有 150 轮 25 次射击的结果，总共 3575 次尝试中有 3750 次命中。

total_shots = 25 * len(results)
total_hits = results.sum()
n, total_shots, total_hits

(25, 3750, np.int64(3575))

因此，整体成功率是 95.3%。

p = total_hits / total_shots
p

np.float64(0.9533333333333334)

现在让我们计算一个Pmf，它代表n=25和刚刚计算出的p值的二项分布。

ps = binomial_pmf(ks, n, p)
pmf_binom = Pmf(ps, ks, name="binomial model")

我们可以将这个分布与实际结果的Pmf进行比较。

pmf_results = Pmf.from_seq(results, name="actual results")

two_bar_plots(pmf_results, pmf_binom)
decorate(xlabel="Hits", ylabel="PMF")

二项模型很好地拟合了数据的分布——尽管它做出了所有竞争者具有相同、不变能力的非现实假设。

泊松分布

作为另一个例子，其中体育赛事的结果遵循可预测的模式，让我们看看冰球比赛中进球的数量。

我们首先模拟一场 60 分钟的比赛，即 3600 秒，假设每场比赛平均进球数为 6 个，进球概率p在任何一秒都是相同的。

n = 3600
m = 6
p = m / 3600
p

0.0016666666666666668

现在我们可以使用以下函数来模拟n秒，并返回进球总数。

def  simulate_goals(n, p):
    return flip(n, p).sum()

如果我们模拟许多场比赛，我们可以确认每场比赛的平均进球数接近 6。

goals = [simulate_goals(n, p) for i in range(1001)]
np.mean(goals)

np.float64(6.021978021978022)

我们可以使用二项分布来模拟这些结果，但当n很大且p很小时，结果也可以很好地由一个泊松分布来模拟，该分布由通常用希腊字母λ表示的值指定，读作lambda，在代码中用变量lam（lambda不是一个合法的变量名，因为它是一个 Python 关键字）表示。lam代表进球率，在示例中为每场比赛 6 个进球。

泊松分布的概率质量函数（PMF）易于计算——给定lam，我们可以使用以下函数来计算一场比赛中出现k个进球的概率。

from  scipy.special  import factorial

def  poisson_pmf(k, lam):
  """Compute the Poisson PMF.

 k (int or array-like): The number of occurrences
 lam (float): The rate parameter (λ) of the Poisson distribution

 returns: float or ndarray
 """
    return (lam**k) * np.exp(-lam) / factorial(k)

SciPy 提供了factorial函数，该函数计算从1到k的整数的乘积。

如果我们用一系列k值调用poisson_pmf，我们可以创建一个Pmf，它表示结果的分布。

lam = 6
ks = np.arange(20)
ps = poisson_pmf(ks, lam)
pmf_poisson = Pmf(ps, ks, name="Poisson model")

并确认分布的均值接近 6。

pmf_poisson.normalize()
pmf_poisson.mean()

np.float64(5.999925498375129)

下图比较了模拟结果与具有相同均值的泊松分布。

pmf_sim = Pmf.from_seq(goals, name="simulation")

two_bar_plots(pmf_sim, pmf_poisson)
decorate(xlabel="Goals", ylabel="PMF")

分布相似，只是由于随机变化存在一些小差异。这并不令人惊讶，因为模拟和泊松模型都是基于相同的假设，即在任何一秒进球的概率是相同的。因此，一个更强的测试是看看模型如何拟合真实数据。

从 HockeyReference 网站上，我下载了 2023-2024 赛季国家冰球联盟（NHL）常规赛中每场比赛的结果（不包括季后赛）。我提取了关于在 60 分钟常规比赛时间内进球的信息，不包括加时赛或点球大战。结果存储在一个 HDF 文件中，每个游戏有一个键，以及一个列表，列出了自比赛开始以来进球的时间，以秒为单位。下载数据的说明在本书的笔记本中。

从 2024 年 7 月 16 日下载的www.hockey-reference.com/leagues/NHL_2024_games.html的原始数据。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/nhl_2023_2024.hdf")

下面是如何从文件中读取键。

filename = "nhl_2023_2024.hdf"

with pd.HDFStore(filename, "r") as store:
    keys = store.keys()

len(keys), keys[0]

(1312, '/202310100PIT')

在常规赛期间共有 1312 场比赛。每个键包含比赛日期和主队的三个字母缩写。我们可以使用read_hdf来查找键并获取进球时间的列表。

times = pd.read_hdf(filename, key=keys[0])
times

0     424
1    1916
2    2137
3    3005
4    3329
5    3513
dtype: int64

在赛季的第一场比赛中，共进了六个球，第一个进球是在比赛进行 424 秒后，最后一个进球是在 3513 秒后——比赛还剩下 87 秒。

3600 - times[5]

np.int64(87)

以下循环读取所有比赛的结果，计算每场比赛的进球数，并将结果存储在列表中。

goals = []

for key in keys:
    times = pd.read_hdf(filename, key=key)
    n = len(times)
    goals.append(n)

每场比赛的平均进球数略超过 6。

lam = np.mean(goals)
lam

np.float64(6.0182926829268295)

我们可以使用poisson_pmf来创建一个Pmf，它代表具有与数据相同均值的泊松分布。

ps = poisson_pmf(ks, lam)
pmf_poisson = Pmf(ps, ks, name="Poisson model")

下面是它与数据的概率质量函数（PMF）相比的图示。

pmf_goals = Pmf.from_seq(goals, name="goals scored")

two_bar_plots(pmf_goals, pmf_poisson)
decorate(xlabel="Goals", ylabel="PMF")

泊松分布很好地拟合了数据，这表明它是冰球进球过程的良好模型。

指数分布

在上一节中，我们模拟了一个简单的冰球比赛模型，其中进球在任何一秒被射中的概率相同。在相同的模型下，结果是，直到第一个进球的时间遵循一个指数分布。

为了演示，让我们再次假设球队平均每场比赛进 6 个球，并计算每秒进球的概率。

n = 3600
m = 6
p = m / 3600
p

0.0016666666666666668

以下函数模拟n秒，并使用argmax找到第一个进球的时间。

def  simulate_first_goal(n, p):
    return flip(n, p).argmax()

这是因为flip的结果是一串 1 和 0，所以最大值几乎总是 1。如果序列中至少有一个进球，argmax返回第一个进球的索引。如果没有进球，它返回 0，但这很少发生，所以我们忽略它。

我们将使用simulate_first_goal模拟 1001 场比赛，并创建一个直到第一个进球的时间列表。

np.random.seed(3)

first_goal_times = [simulate_first_goal(n, p) for i in range(1001)]
mean = np.mean(first_goal_times)
mean

np.float64(597.7902097902098)

直到第一个进球的平均时间是接近 600 秒，即 10 分钟。这很有道理——如果我们预计每 60 分钟比赛有 6 个进球，那么平均每 10 分钟就会有一个进球。

当n很大且p很小时，我们可以从数学上证明第一个进球的期望时间遵循指数分布。

因为模拟生成了许多独特的时间值，我们将使用累积分布函数（CDFs）来比较分布，而不是概率质量函数（PMFs）。指数分布的 CDF 很容易计算。

def  exponential_cdf(x, lam):
  """Compute the exponential CDF.

 x: float or sequence of floats
 lam: rate parameter

 returns: float or NumPy array of cumulative probability
 """
    return 1 - np.exp(-lam * x)

lam的值是单位时间内的平均事件数——在这个例子中是每秒的进球数。我们可以使用模拟结果的平均值来计算lam。

lam = 1 / mean
lam

np.float64(0.0016728276636563566)

如果我们用时间值的范围调用此函数，我们可以近似第一次进球时间的分布。NumPy 函数linspace创建一个等间距值的数组；在这个例子中，它计算从 0 到 3600 的 201 个值，包括两者。

from  empiricaldist  import Cdf

ts = np.linspace(0, 3600, 201)
ps = exponential_cdf(ts, lam)
cdf_expo = Cdf(ps, ts, name="exponential model")

下图比较了模拟结果与我们刚刚计算的指数分布。

cdf_sim = Cdf.from_seq(first_goal_times, name="simulation")

cdf_expo.plot(ls=":", color="gray")
cdf_sim.plot()

decorate(xlabel="Time of first goal (seconds)", ylabel="CDF")

指数模型与模拟结果非常吻合——但更严格的测试是看看它在真实数据上的表现如何。

下面的循环读取所有比赛的结果，获取第一次进球的时间，并将结果存储在列表中。如果没有进球，它会在列表中添加nan。

filename = "nhl_2023_2024.hdf"

with pd.HDFStore(filename, "r") as store:
    keys = store.keys()

firsts = []

for key in keys:
    times = pd.read_hdf(filename, key=key)
    if len(times) > 0:
        firsts.append(times[0])
    else:
        firsts.append(np.nan)

为了估计进球率，我们可以使用nanmean，它计算时间的平均值，忽略nan值。

lam = 1 / np.nanmean(firsts)
lam

np.float64(0.0015121567467720825)

现在，我们可以使用与数据相同的进球率来计算指数分布的累积分布函数（CDF）。

ps = exponential_cdf(ts, lam)
cdf_expo = Cdf(ps, ts, name="exponential model")

为了计算数据的 CDF，我们将使用dropna=False参数，它包括末尾的nan值。

cdf_firsts = Cdf.from_seq(firsts, name="data", dropna=False)
cdf_firsts.tail()

	probs
3286.0	0.996951
3581.0	0.997713
NaN	1.000000

下图比较了指数分布与数据的分布。

cdf_expo.plot(ls=":", color="gray")
cdf_firsts.plot()

decorate(xlabel="Time of first goal (seconds)", ylabel="CDF")

数据在某些地方偏离了模型——看起来在模型预测的前 1000 秒内进球比模型预测的少。但仍然，模型很好地拟合了数据。

这些模型的基础假设——进球的泊松模型和时间指数模型——是进球在任何游戏秒数中发生的可能性是相同的。如果你问一个曲棍球迷这是否正确，他们会说不是，而且他们会是对的——现实世界在许多方面违反了这样的假设。尽管如此，理论分布通常与真实数据非常吻合。

正态分布

我们在现实世界中测量的许多事物都遵循正态分布，也称为高斯分布或“钟形曲线”。为了了解这些分布的来源，让我们考虑一个关于巨型南瓜生长的模型。假设每天，如果天气不好，南瓜会增重 1 磅；如果天气晴朗，增重 2 磅；如果天气好，增重 3 磅。并且假设每天的天气以相同的概率是不好、晴朗或好。

我们可以使用以下函数来模拟n天的模型并返回总重量增加量。

def  simulate_growth(n):
    choices = [1, 2, 3]
    gains = np.random.choice(choices, n)
    return gains.sum()

NumPy 的random模块提供了一个choice函数，该函数从一个值序列中生成n个随机选择，在这个例子中是choices。

现在假设有 1001 人在不同地方、不同天气下种植巨型南瓜。如果我们模拟 100 天的生长过程，我们得到一个包含 1001 个重量的列表。

sim_weights = [simulate_growth(100) for i in range(1001)]
m, s = np.mean(sim_weights), np.std(sim_weights)
m, s

(np.float64(199.37062937062936), np.float64(8.388630840376777))

平均值接近 200 磅，标准差约为 8 磅。为了查看重量是否遵循正态分布，我们将使用以下函数，该函数接受一个样本并创建一个Cdf，它代表具有与样本相同均值和标准差的正态分布，评估范围从平均值以下 4 个标准差到平均值以上 4 个标准差。

from  scipy.stats  import norm

def  make_normal_model(data):
    m, s = np.mean(data), np.std(data)
    low, high = m - 4 * s, m + 4 * s
    qs = np.linspace(low, high, 201)
    ps = norm.cdf(qs, m, s)
    return Cdf(ps, qs, name="normal model")

这就是我们的使用方法。

cdf_model = make_normal_model(sim_weights)

现在，我们可以创建一个代表模拟结果分布的Cdf。

cdf_sim_weights = Cdf.from_seq(sim_weights, name="simulation")

我们将使用以下函数来比较分布。cdf_model和cdf_data是Cdf对象。xlabel是一个字符串，options是一个选项字典，它控制cdf_data的绘图方式。

def  two_cdf_plots(cdf_model, cdf_data, xlabel="", **options):
    cdf_model.plot(ls=":", color="gray")
    cdf_data.plot(**options)
    decorate(xlabel=xlabel, ylabel="CDF")

下面是结果。

two_cdf_plots(cdf_model, cdf_sim_weights, xlabel="Weight (pounds)")

正态模型很好地拟合了权重的分布。一般来说，当我们累加足够的随机因素时，总和往往会遵循正态分布。这是中心极限定理的结果，我们将在第十四章中再次讨论。

但首先让我们看看正态分布如何拟合实际数据。作为一个例子，我们将查看国家家庭增长调查（NSFG）中出生体重的分布。我们可以使用read_fem_preg读取数据，然后选择totalwgt_lb列，该列记录以磅为单位的出生体重。

以下单元格下载数据文件并安装statadict，我们需要它来读取数据。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz")

try:
    import  statadict
except ImportError:
    %pip install statadict

import  nsfg

preg = nsfg.read_fem_preg()
birth_weights = preg["totalwgt_lb"].dropna()

出生体重的平均值约为 7.27 磅，标准差为 1.4 磅，但正如我们所看到的，这个数据集中有一些异常值，可能是错误。

m, s = np.mean(birth_weights), np.std(birth_weights)
m, s

(np.float64(7.265628457623368), np.float64(1.40821553384062))

为了减少异常值对估计均值和标准差的影响，我们将使用 SciPy 函数trimboth来移除最高和最低值。

from  scipy.stats  import trimboth

trimmed = trimboth(birth_weights, 0.01)
m, s = np.mean(trimmed), np.std(trimmed)
m, s

(np.float64(7.280883100022579), np.float64(1.2430657948614345))

使用修剪后的数据，平均值略低，标准差显著降低。我们将使用修剪后的数据来制作正态模型。

cdf_model = make_normal_model(trimmed)

并将其与数据的Cdf进行比较。

cdf_birth_weight = Cdf.from_seq(birth_weights, name='data')

two_cdf_plots(cdf_model, cdf_birth_weight, xlabel="Birth weight (pounds)")

正态模型很好地拟合了数据，除了 5 磅以下的部分，那里的数据分布位于模型左侧——也就是说，最轻的婴儿比我们预期的正态分布要轻。现实世界通常比简单的数学模型更复杂。

在上一节中，我们假设南瓜每天增长 1-3 磅，根据天气情况模拟南瓜的生长。相反，让我们假设它们的增长与当前重量成正比，因此大南瓜每天增长的重量比小南瓜多——这可能是更现实的。

以下函数模拟这种比例增长，其中南瓜在天气恶劣时增加其体重的 3%，在天气晴朗时增加 5%，在天气良好时增加 7%。同样，我们将假设在任何给定的一天，天气恶劣、晴朗或良好的概率是相等的。

def  simulate_proportionate_growth(n):
    choices = [1.03, 1.05, 1.07]
    gains = np.random.choice(choices, n)
    return gains.prod()

如果南瓜增加了 3%的体重，最终体重是初始体重与因子 1.03 的乘积。因此，我们可以通过选择随机因子并将它们相乘来计算 100 天后的体重。

我们将调用这个函数 1001 次来模拟 1001 个南瓜并保存它们的体重。

sim_weights = [simulate_proportionate_growth(100) for i in range(1001)]
np.mean(sim_weights), np.std(sim_weights)

(np.float64(130.80183363824722), np.float64(20.956047434921466))

平均体重约为 131 磅；标准差约为 21 磅。因此，在这个模型中，南瓜比上一个模型小，但变化更大。

我们可以用数学方法证明它们遵循对数正态分布，这意味着体重的对数遵循正态分布。为了检查，我们将计算体重的对数以及它们的均值和标准差。我们可以使用任何底数的对数，但我将使用底数为 10，因为这使结果更容易解释。

log_sim_weights = np.log10(sim_weights)
m, s = np.mean(log_sim_weights), np.std(log_sim_weights)
m, s

(np.float64(2.1111299372609933), np.float64(0.06898607064749827))

现在我们将比较对数分布与具有相同均值和标准差的正态分布。

cdf_model = make_normal_model(log_sim_weights)
cdf_log_sim_weights = Cdf.from_seq(log_sim_weights, name="simulation")

two_cdf_plots(
    cdf_model, cdf_log_sim_weights, xlabel="Pumpkin weight (log10 pounds)"
)

模型与模拟结果非常吻合，这正是我们预期的。

如果人们的体重变化与南瓜相似，即每年的体重变化与当前体重成比例，我们可能会预期成年体重的分布遵循对数正态分布。让我们找出答案。

国家慢性疾病预防和健康促进中心作为行为风险因素监测系统（BRFSS）的一部分进行年度调查。2008 年，他们采访了 414,509 名受访者，并询问了他们的人口统计信息、健康状况和健康风险。他们收集的数据中包括 398,484 名受访者的体重。下载数据的说明在本书的笔记本中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/CDBRFS08.ASC.gz")

thinkstats模块提供了一个函数，用于读取 BRFSS 数据并返回 Pandas DataFrame。

from  thinkstats  import read_brfss

brfss = read_brfss()

成人体重以公斤为单位记录在wtkg2列中。

adult_weights = brfss["wtkg2"].dropna()
m, s = np.mean(adult_weights), np.std(adult_weights)
m, s

(np.float64(78.9924529968581), np.float64(19.546132387397257))

平均体重约为 79 公斤。在我们计算对数之前，让我们看看体重是否遵循正态分布。

cdf_model = make_normal_model(adult_weights)
cdf_adult_weights = Cdf.from_seq(adult_weights, name="adult weight")

two_cdf_plots(cdf_model, cdf_adult_weights, xlabel="Adult weight (kilograms)")

对于某些目的来说，正态分布可能足够好地描述这些数据，但让我们看看我们是否能做得更好。

这是对数转换后体重的分布和具有相同均值和标准差的正态模型。

log_adult_weights = np.log10(adult_weights)
cdf_model = make_normal_model(log_adult_weights)

cdf_log_adult_weights = Cdf.from_seq(log_adult_weights, name="log adult weight")

two_cdf_plots(cdf_model, cdf_log_adult_weights, xlabel="Adult weight (log10 kg)")

正态模型对对数的拟合优于对体重的拟合，这表明比例增长是体重增加比加性增长更好的模型。

为什么建模？

在本章的开头，我说过许多现实世界现象可以用理论分布来模拟。但可能不清楚为什么我们应该关心这一点。

与所有模型一样，理论分布是抽象的，这意味着它们忽略了被认为无关的细节。例如，观察到的分布可能包含测量误差或特定于样本的怪癖；理论模型忽略了这些特殊性。

理论模型也是一种数据压缩形式。当一个模型很好地拟合数据集时，一组小的数字可以总结大量数据。

当自然现象的数据符合理论分布时，有时会令人惊讶，但这些观察可以提供对物理系统的洞察。有时我们可以解释为什么观察到的分布具有特定的形式。例如，在前一节中我们发现成年人的体重很好地符合对数正态分布，这表明体重年复一年的变化可能与当前体重成比例。

此外，理论分布适合进行数学分析，正如我们在第十四章 Chapter 14 中将会看到的。

但重要的是要记住，所有模型都是不完美的。现实世界的数据永远不会完美地符合理论分布。人们有时会像谈论数据是由模型生成的；例如，他们可能会说人类身高的分布是正态的，或收入的分布是对数正态的。字面意义上，这些说法不能成立——现实世界和数学模型之间总是存在差异。

如果模型能够捕捉现实世界的相关方面并省略不必要的细节，那么模型是有用的。但相关或不必要取决于你打算如何使用该模型。

术语表

二项分布：一种理论分布，常用于模拟一系列成功或击中的次数。
泊松分布：一种理论分布，常用于模拟在一段时间内发生的事件数量。
指数分布：一种理论分布，常用于模拟事件之间的时间间隔。
正态分布：一种理论分布，常用于模拟呈对称、钟形曲线的数据。
对数正态分布：一种理论分布，常用于模拟遵循向右偏斜的钟形曲线的数据。

练习

练习 5.1

在 NSFG 受访者文件中，numfmhh列记录了每个受访者的家庭“家庭成员数量”。我们可以使用read_fem_resp来读取文件，并使用query来选择在受访时年龄为 25 岁或以上的受访者。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dat.gz")

from  nsfg  import read_fem_resp

resp = read_fem_resp()

older = resp.query("age >= 25")
num_family = older["numfmhh"]

计算这些较老受访者的 numfmhh 的概率质量函数（Pmf），并将其与具有相同均值的泊松分布进行比较。泊松模型与数据的拟合程度如何？

练习 5.2

在本章的早期，我们看到了在曲棍球比赛中第一个进球的时间遵循指数分布。如果我们的进球模型是正确的，那么在任何时间进球的可能性都是相同的，无论上一个进球已经过去多久。如果这是真的，我们预计进球之间的时间间隔也将遵循指数分布。

以下循环再次读取曲棍球数据，计算连续进球之间的时间（如果一场比赛中有多于一个进球），并将进球间时间收集到一个列表中。

filename = "nhl_2023_2024.hdf"

with pd.HDFStore(filename, "r") as store:
    keys = store.keys()

intervals = []

for key in keys:
    times = pd.read_hdf(filename, key=key)
    if len(times) > 1:
        intervals.extend(times.diff().dropna())

使用 exponential_cdf 计算与观察到的间隔具有相同均值的指数分布的累积分布函数（CDF），并将此模型与数据的 CDF 进行比较。

练习 5.3

人类身高的分布更像是正态分布还是对数正态分布？为了找出答案，我们可以从 BRFSS 中选择身高数据，如下所示：

adult_heights = brfss["htm3"].dropna()
m, s = np.mean(adult_heights), np.std(adult_heights)
m, s

(np.float64(168.82518961012298), np.float64(10.35264015645592))

计算这些值的累积分布函数（CDF），并将其与具有相同均值和标准差的正态分布进行比较。然后计算高度的对数，并将对数分布与正态分布进行比较。根据视觉比较，哪个模型更适合数据？

《Think Stats：Python 中的探索性数据分析》（第 3 版）

版权 2024 艾伦·B·唐尼

代码许可：MIT License

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

概率密度函数

原文：allendowney.github.io/ThinkStats/chap06.html

在上一章中，我们使用二项分布、泊松分布、指数分布和正态分布等理论分布来模拟数据。

二项分布和泊松分布是离散的，这意味着结果必须是不同的或独立的元素，比如击中和失误的整数次数，或者进球数。在离散分布中，每个结果都与一个概率质量相关联。

指数分布和正态分布是连续的，这意味着结果可以位于可能值范围内的任何点。在连续分布中，每个结果都与一个概率密度相关联。概率密度是一个抽象的概念，许多人一开始会觉得很难理解，但我们会一步一步来。作为第一步，让我们再次思考如何比较分布。

点击此处运行此笔记本在 Colab 上。

from  os.path  import basename, exists

def  download(url):
    filename = basename(url)
    if not exists(filename):
        from  urllib.request  import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("Downloaded " + local)

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")

try:
    import  empiricaldist
except ImportError:
    %pip install empiricaldist

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  plt

from  thinkstats  import decorate

比较分布

在上一章中，当我们比较离散分布时，我们使用条形图来展示它们的概率质量函数（PMFs）。当我们比较连续分布时，我们使用线形图来展示它们的累积分布函数（CDFs）。

对于离散分布，我们也可以使用累积分布函数（CDFs）。例如，这里有一个参数为lam=2.2的泊松分布的概率质量函数（PMF），这是一个很好的模型，用于描述 NSFG 数据中家庭规模的分布。

下面的单元格下载数据文件并安装statadict，这是我们读取数据所需的。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/nsfg.py")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemResp.dat.gz")

try:
    import  statadict
except ImportError:
    %pip install statadict

我们可以使用read_fem_resp来读取受访者数据文件。

from  nsfg  import read_fem_resp

resp = read_fem_resp()

接下来，我们将选择 25 岁及以上人群的家庭规模。

older = resp.query("age >= 25")
num_family = older["numfmhh"]

然后创建一个Pmf来表示响应的分布。

from  empiricaldist  import Pmf

pmf_family = Pmf.from_seq(num_family, name="data")

这里还有一个表示具有相同均值的泊松分布的Pmf。

from  thinkstats  import poisson_pmf

lam = 2.2
ks = np.arange(11)
ps = poisson_pmf(ks, lam)

pmf_poisson = Pmf(ps, ks, name="Poisson model")

下面是数据分布与泊松模型的比较。

from  thinkstats  import two_bar_plots

two_bar_plots(pmf_family, pmf_poisson)
decorate(xlabel="Number of family members")

比较 PMFs，我们可以看到模型很好地拟合了数据，但有一些偏差。

为了了解这些偏差的实质，比较 CDFs 可能会有所帮助。我们可以使用make_cdf来计算数据和模型的 CDFs。

cdf_family = pmf_family.make_cdf()
cdf_poisson = pmf_poisson.make_cdf()

下面是它们的形状。

from  thinkstats  import two_cdf_plots

two_cdf_plots(cdf_poisson, cdf_family)
decorate(xlabel="Number of family members")

当我们比较 CDFs 时，偏差不那么明显，但我们能看到分布差异的具体位置和方式。概率质量函数（PMFs）往往强调小的差异——有时 CDFs 能更好地展现整体情况。

CDFs 也适用于连续数据。作为一个例子，让我们再次看看出生体重的分布，它位于 NSFG 怀孕文件中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dct")
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/2002FemPreg.dat.gz")

from  nsfg  import read_fem_preg

preg = read_fem_preg()
birth_weights = preg["totalwgt_lb"].dropna()

这是我们在上一章中使用来拟合数据的正态模型的代码。

from  scipy.stats  import trimboth
from  thinkstats  import make_normal_model

trimmed = trimboth(birth_weights, 0.01)
cdf_model = make_normal_model(trimmed)

下面是数据分布与正态模型的比较。

from  empiricaldist  import Cdf

cdf_birth_weight = Cdf.from_seq(birth_weights, name="sample")
two_cdf_plots(cdf_model, cdf_birth_weight, xlabel="Birth weight (pounds)")

正如我们在上一章中看到的，正态模型很好地拟合了数据，除了最轻婴儿的范围。

在我看来，累积分布函数（CDFs）通常是比较数据与模型的最佳方式。但对于不熟悉 CDFs 的观众来说，还有一个选项：概率密度函数。

概率密度

我们将从正态分布的概率密度函数（PDF）开始，该函数计算给定mu和sigma的量xs的密度。

def  normal_pdf(xs, mu, sigma):
  """Evaluates the normal probability density function."""
    z = (xs - mu) / sigma
    return np.exp(-(z**2) / 2) / sigma / np.sqrt(2 * np.pi)

对于mu和sigma，我们将使用修剪后新生儿体重的均值和标准差。

m, s = np.mean(trimmed), np.std(trimmed)

现在我们将评估一系列重量的normal_pdf。

low = m - 4 * s
high = m + 4 * s
qs = np.linspace(low, high, 201)
ps = normal_pdf(qs, m, s)

然后绘制它。

plt.plot(qs, ps, label="normal model", ls=":", color="gray")
decorate(xlabel="Birth weight (pounds)", ylabel="Density")

结果看起来像钟形曲线，这是正态分布的特征。

当我们评估normal_pdf时，结果是概率密度。例如，这里是在均值处评估的密度函数，这是密度最高的地方。

normal_pdf(m, m, s)

np.float64(0.32093416297880123)

单独来看，概率密度并没有什么意义——最重要的是，它不是一个概率。说随机选择的新生儿体重等于m的概率是 32%是不正确的。事实上，新生儿体重真正、完全、精确地等于m——或任何其他特定值——的概率是零。

然而，我们可以使用概率密度来计算结果落在两个值之间的区间的概率，通过计算曲线下的面积。

我们可以用normal_pdf函数做到这一点，但使用定义在thinkstats模块中的NormalPdf类更方便。以下是如何创建与 NSFG 数据集中新生儿体重相同的均值和标准差的NormalPdf对象。

from  thinkstats  import NormalPdf

pdf_model = NormalPdf(m, s, name="normal model")
pdf_model

NormalPdf(7.280883100022579, 1.2430657948614345, name='normal model')

如果我们像调用函数一样调用这个对象，它将评估正态 PDF。

pdf_model(m)

np.float64(0.32093416297880123)

现在，为了计算概率密度函数（PDF）下的面积，我们可以使用以下函数，该函数接受一个NormalPdf对象和一个区间的边界，low和high。它评估low和high之间等间距的量，并使用 SciPy 函数simpson来估计曲线下的面积（simpson之所以这样命名，是因为它使用了一个称为辛普森方法的算法）。

from  scipy.integrate  import simpson

def  area_under(pdf, low, high):
    qs = np.linspace(low, high, 501)
    ps = pdf(qs)
    return simpson(y=ps, x=qs)

如果我们从图中的最低点到最高点计算曲线下的面积，结果接近 1。

area_under(pdf_model, 2, 12)

np.float64(0.9999158086616793)

如果我们将区间从负无穷大到正无穷大扩展，总面积正好是 1。

如果我们从 0——或任何远低于均值的值开始——我们可以计算体重小于或等于 8.5 磅的婴儿体重的比例。

area_under(pdf_model, 0, 8.5)

np.float64(0.8366380335513807)

你可能还记得，“小于或等于给定值的分数”是累积分布函数（CDF）的定义。因此，我们可以使用正态分布的 CDF 来计算相同的结果。

from  scipy.stats  import norm

norm.cdf(8.5, m, s)

np.float64(0.8366380358092718)

同样，我们可以使用密度曲线下的面积来计算 6 到 8 磅之间出生体重的分数。

area_under(pdf_model, 6, 8)

np.float64(0.5671317752927691)

或者，我们可以使用 CDF 来计算小于 8 的分数，然后减去小于 6 的分数。

norm.cdf(8, m, s) - norm.cdf(6, m, s)

np.float64(0.5671317752921801)

因此，CDF 是 PDF 曲线下的面积。如果你知道微积分，另一种表达相同意思的方式是 CDF 是 PDF 的积分。反之，PDF 是 CDF 的导数。

指数概率密度函数（PDF）

为了理解概率密度，看看另一个例子可能会有所帮助。在前一章中，我们使用指数分布来模拟冰球比赛中第一个进球的时间。我们使用了以下函数来计算指数 CDF，其中lam是每单位时间进球的速率。

def  exponential_cdf(x, lam):
  """Compute the exponential CDF.

 x: float or sequence of floats
 lam: rate parameter

 returns: float or NumPy array of cumulative probability
 """
    return 1 - np.exp(-lam * x)

我们可以这样计算指数分布的 PDF。

def  exponential_pdf(x, lam):
  """Evaluates the exponential PDF.

 x: float or sequence of floats
 lam: rate parameter

 returns: float or NumPy array of probability density
 """
    return lam * np.exp(-lam * x)

thinkstats提供了一个使用此函数来计算指数 PDF 的ExponentialPdf对象。我们可以用它来表示每场比赛 6 个进球的指数分布。

from  thinkstats  import ExponentialPdf

lam = 6
pdf_expo = ExponentialPdf(lam, name="exponential model")
pdf_expo

ExponentialPdf(6, name='exponential model')

ExponentialPdf提供了一个plot函数，我们可以用它来绘制 PDF - 注意，这里的单位是比赛，而不是前一章中的秒。

qs = np.linspace(0, 1.5, 201)
pdf_expo.plot(qs, ls=":", color="gray")
decorate(xlabel="Time (games)", ylabel="Density")

观察 y 轴，你可能注意到这些密度中的一些大于 1，这是一个提醒，概率密度不是一个概率。但是，密度曲线下的面积是一个概率，因此它永远不会大于 1。

如果我们从 0 到 1.5 场比赛的曲线下计算面积，我们可以确认结果接近 1。

area_under(pdf_expo, 0, 1.5)

np.float64(0.999876590779019)

如果我们将区间扩展得更大，结果略大于 1，但这是因为我们在数值上近似面积。从数学上讲，它正好是 1，正如我们可以使用指数 CDF 来确认的那样。

from  thinkstats  import exponential_cdf

exponential_cdf(7, lam)

np.float64(1.0)

我们可以使用密度曲线下的面积来计算任何区间内进球的概率。例如，这是 60 分钟比赛中第一分钟进球的概率。

area_under(pdf_expo, 0, 1 / 60)

np.float64(0.09516258196404043)

我们可以使用指数 CDF 来计算相同的结果。

exponential_cdf(1 / 60, lam)

np.float64(0.09516258196404048)

总结来说，如果我们评估 PDF，结果是概率密度 - 这不是一个概率。然而，如果我们计算 PDF 下的面积，结果是某个量落在区间内的概率。或者，我们也可以通过在区间的开始和结束处评估 CDF 并计算差值来找到相同的概率。

比较概率质量函数（PMFs）和概率密度函数（PDFs）

将样本的 PMF 与理论模型的 PDF 进行比较是一个常见的错误。例如，假设我们想将出生体重的分布与正态模型进行比较。这是一个表示数据分布的Pmf。

pmf_birth_weight = Pmf.from_seq(birth_weights, name="data")

我们已经有了pdf_model，它代表具有相同均值和标准差的正态分布的 PDF。如果我们把它们画在同一个轴上，会发生什么。

pdf_model.plot(ls=":", color="gray")
pmf_birth_weight.plot()

decorate(xlabel="Birth weight (pounds)", ylabel="PMF? Density?")

它效果不太好。一个原因是它们不在同一个单位上。PMF 包含概率质量，而 PDF 包含概率密度，所以我们不能比较它们，也不应该在同一个轴上绘制它们。

作为解决这个问题的第一次尝试，我们可以制作一个Pmf，通过在离散的点集上评估 PDF 来近似正态分布。NormalPdf提供了一个make_pmf方法来做这件事。

pmf_model = pdf_model.make_pmf()

结果是一个包含概率质量的标准化Pmf，所以我们至少可以把它画在和数据 PMF 相同的轴上。

pmf_model.plot(ls=":", color="gray")
pmf_birth_weight.plot()

decorate(xlabel="Birth weight (pounds)", ylabel="PMF")

但这仍然不是比较分布的好方法。一个问题就是两个Pmf对象包含不同数量的量，而且pmf_birth_weight中的量不是等间距的，所以概率质量实际上是不可比较的。

len(pmf_model), len(pmf_birth_weight)

(201, 184)

另一个问题就是数据的Pmf是嘈杂的。所以让我们试试别的。

核密度估计

我们不是用模型来制作 PMF，而是可以用数据来制作 PDF。为了展示这是如何工作的，我将从一个小的数据样本开始。

# Set the random seed so we get the same results every time
np.random.seed(3)

n = 10
sample = birth_weights.sample(n)

这个样本的Pmf看起来是这样的。

for weight in sample:
    pmf = Pmf.from_seq([weight]) / n
    pmf.bar(width=0.08, alpha=0.5)

xlim = [1.5, 12.5]
decorate(xlabel="Birth weight (pounds)", ylabel="PMF", xlim=xlim)

这种表示分布的方式将数据视为离散的，所以每个概率质量都堆叠在单个点上。但出生重量实际上是连续的，所以测量之间的量也是可能的。我们可以通过用连续概率密度替换每个离散概率质量来表示这种可能性，如下所示。

qs = np.linspace(2, 12, 201)

for weight in sample:
    ps = NormalPdf(weight, 0.75)(qs) / n
    plt.plot(qs, ps, alpha=0.5)

decorate(xlabel="Birth weight (pounds)", ylabel="PDF", xlim=xlim)

对于样本中的每个重量，我们创建一个以观察到的重量为均值的NormalPdf – 现在我们来把它们加起来。

low_ps = np.zeros_like(qs)

for weight in sample:
    ps = NormalPdf(weight, 0.75)(qs) / n
    high_ps = low_ps + ps
    plt.fill_between(qs, low_ps, high_ps, alpha=0.5, lw=1, ec="white")
    low_ps = high_ps

decorate(xlabel="Birth weight (pounds)", ylabel="PDF", xlim=xlim)

当我们为每个数据点加总概率密度时，结果是整个样本的概率密度估计。这个过程被称为核密度估计或 KDE。在这个上下文中，“核”是我们加总的小密度函数之一。因为我们使用的核是正态分布，也称为高斯分布，所以我们更具体地说，我们计算了一个高斯 KDE。

SciPy 提供了一个名为 gaussian_kde 的函数，该函数实现了此算法。以下是我们可以如何使用它来估计出生体重的分布。

from  scipy.stats  import gaussian_kde

kde = gaussian_kde(birth_weights)

结果是一个表示估计 PDF 的对象，我们可以通过像调用函数一样调用它来评估它。

ps = kde(qs)

这就是结果的样子。

plt.plot(qs, ps)

decorate(xlabel="Birth weight (pounds)", ylabel="Density")

thinkstats 提供了一个 Pdf 对象，它接受 gaussian_kde 的结果，以及一个表示应在何处评估密度的域。以下是我们的实现方法。

from  thinkstats  import Pdf

domain = np.min(birth_weights), np.max(birth_weights)
kde_birth_weights = Pdf(kde, domain, name="data")

Pdf 提供了一个 plot 方法，我们可以用它来比较样本的估计 PDF 与正态分布的 PDF。

pdf_model.plot(ls=":", color="gray")
kde_birth_weights.plot()

decorate(xlabel="Birth weight (pounds)", ylabel="Density")

核密度估计使得将数据集的分布与理论模型进行比较成为可能，对于某些听众来说，这是一种很好的可视化比较的方法。但对于熟悉 CDFs 的听众来说，比较 CDFs 通常更好。

分布框架

到目前为止，我们有一套完整的方式来表示分布：PMFs、CDFs 和 PDFs。以下图显示了这些表示以及从一个到另一个的转换。例如，如果我们有一个 Pmf，我们可以使用 cumsum 函数来计算概率的累积和，从而得到一个表示相同分布的 Cdf。

链接

为了演示这些转换，我们将使用一个新的数据集，该数据集“包含在澳大利亚布里斯班一家医院 24 小时内出生的 44 个婴儿的出生时间、性别和出生体重”，根据描述。下载数据的说明在本书的笔记本中。

根据文件中的信息

来源：Steele, S. (1997 年 12 月 21 日)，“圣诞节的十二个婴儿：24 小时婴儿潮”，The Sunday Mail（布里斯班），第 7 页。

数据背后的故事：1997 年 12 月 18 日，在澳大利亚昆士兰州布里斯班 Mater Mothers’ Hospital，在 24 小时内出生了 44 个婴儿——创下了新纪录。对于这 44 个婴儿中的每一个，The Sunday Mail 记录了出生时间、孩子的性别和出生体重（克）。

关于此数据集的更多信息可以在《统计学教育杂志》中找到“数据集和故事”文章“用于演示常见分布的简单数据集”（Dunn 1999）。

从jse.amstat.org/datasets/babyboom.txt下载

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/babyboom.dat")

我们可以像这样读取数据。

from  thinkstats  import read_baby_boom

boom = read_baby_boom()
boom.head()

	时间	性别	体重 _g	分钟
0	5	1	3837	5
1	104	1	3334	64
2	118	2	3554	78
3	155	2	3838	115
4	257	2	3625	177

minutes 列记录“每个出生自午夜以来的分钟数”。因此，我们可以使用 diff 方法来计算每个连续出生之间的间隔。

diffs = boom["minutes"].diff().dropna()

如果出生在任何一分钟内发生的概率相等，我们预计这些间隔将遵循指数分布。在现实中，这个假设并不完全准确，但指数分布可能仍然是一个好的数据模型。

为了找出答案，我们将首先创建一个表示间隔分布的 Pmf。

pmf_diffs = Pmf.from_seq(diffs, name="data")
pmf_diffs.bar(width=1)

decorate(xlabel="Interval (minutes)", ylabel="PMF")

然后，我们可以使用 make_cdf 计算累积概率并将它们存储在 Cdf 对象中。

cdf_diffs = pmf_diffs.make_cdf()
cdf_diffs.step()

decorate(xlabel="Interval (minutes)", ylabel="CDF")

Pmf 和 Cdf 在意义上是等效的，即如果我们给出任何一个，我们都可以计算出另一个。为了演示，我们将使用 make_pmf 方法，该方法计算 Cdf 中连续概率之间的差异，并返回一个 Pmf。

pmf_diffs2 = cdf_diffs.make_pmf()

结果应该与原始 Pmf 相同，但可能会有小的浮点数错误。我们可以使用 allclose 来检查结果是否接近原始 Pmf。

np.allclose(pmf_diffs, pmf_diffs2)

True

它确实是。

从一个 Pmf 中，我们可以通过使用 Pmf 中的概率作为权重调用 gaussian_kde 来估计一个密度函数。

kde = gaussian_kde(pmf_diffs.qs, weights=pmf_diffs.ps)

要绘制结果，我们可以使用 kde 创建一个 Pdf 对象，并调用 plot 方法。

domain = np.min(pmf_diffs.qs), np.max(pmf_diffs.qs)
kde_diffs = Pdf(kde, domain=domain, name="estimated density")

kde_diffs.plot(ls=":", color="gray")
decorate(xlabel="Interval (minutes)", ylabel="Density")

为了查看估计的密度是否遵循指数模型，我们可以使用与数据相同均值的 ExponentialCdf。

from  thinkstats  import ExponentialCdf

m = diffs.mean()
lam = 1 / m
cdf_model = ExponentialCdf(lam, name="exponential CDF")

这是与数据 CDF 相比的样子。

cdf_model.plot(ls=":", color="gray")
cdf_diffs.step()

decorate(xlabel="Interval (minutes)", ylabel="CDF")

指数模型很好地拟合了数据的 CDF。

给定一个 ExponentialCdf，我们可以使用 make_cdf 来离散化CDF – 即通过在一系列等间距的数量上评估 CDF 来得到一个离散近似。

qs = np.linspace(0, 160)
discrete_cdf_model = cdf_model.make_cdf(qs)
discrete_cdf_model.step(color="gray")

decorate(xlabel="Interval (minutes)", ylabel="CDF")

最后，要从离散 CDF 转换到连续 CDF，我们可以在步骤之间进行插值，这就是如果我们使用 plot 方法而不是 step 方法时看到的情况。

discrete_cdf_model.plot(color="gray")

decorate(xlabel="Interval (minutes)", ylabel="CDF")

最后，PDF 是连续 CDF 的导数，而 CDF 是 PDF 的积分。

为了演示，我们可以使用 SymPy 定义指数分布的累积分布函数（CDF）并计算其导数。

import  sympy  as  sp

x = sp.Symbol("x", real=True, positive=True)
λ = sp.Symbol("λ", real=True, positive=True)

cdf = 1 - sp.exp(-λ * x)
cdf

[\displaystyle 1 - e^{- x λ}]

pdf = sp.diff(cdf, x)
pdf

[\displaystyle λ e^{- x λ}]

如果我们对结果进行积分，我们就会得到 CDF，尽管在这个过程中我们失去了积分常数。

sp.integrate(pdf, x)

[\displaystyle - e^{- x λ}]

这个例子展示了我们如何使用 Pmf、Cdf 和 Pdf 对象来表示概率质量函数、累积分布函数和概率密度函数，并演示了将每个转换为其他的过程。

术语表

连续的：如果一个量可以在数轴上的一个范围内取任何值，则该量是连续的。我们世界中测量的大多数事物——如重量、距离和时间——都是连续的。
离散的：如果一个量可以有一个有限的值集，如整数或类别，则该量是离散的。精确计数和分类变量也是离散的。
概率密度函数 (PDF)：一个函数，它显示了密度（不是概率）是如何在连续变量的值上分布的。PDF 在区间下的面积给出了变量落在该区间范围内的概率。
概率密度：PDF 在特定点的值；它本身不是概率，但它可以用来计算概率。
核密度估计 (KDE)：一种基于样本估计 PDF 的方法。
离散化：通过将连续量的范围划分为离散的级别或类别来近似连续量。

练习

练习 6.1

在世界杯足球（足球）比赛中，假设直到第一个进球的时间可以用速率 lam=2.5 个进球每场比赛的指数分布很好地建模。创建一个 ExponentialPdf 来表示这个分布，并使用 area_under 来计算直到第一个进球的时间小于半场比赛时间的概率。然后使用 ExponentialCdf 来计算相同的概率，并检查结果是否一致。

使用 ExponentialPdf 来计算第一个进球在比赛下半场得分的概率。然后使用 ExponentialCdf 来计算相同的概率，并检查结果是否一致。

练习 6.2

为了加入蓝人乐队，你必须是一名身高在 5’10” 到 6’1” 之间的男性，这大约是 178 到 185 厘米。让我们看看在美国的成年男性人口中有多少比例符合这一要求。

BRFSS 中男性参与者的身高可以用均值为 178 厘米和标准差 7 厘米的正态分布很好地建模。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/CDBRFS08.ASC.gz")

from  thinkstats  import read_brfss

brfss = read_brfss()
male = brfss.query("sex == 1")
heights = male["htm3"].dropna()

from  scipy.stats  import trimboth

trimmed = trimboth(heights, 0.01)
m, s = np.mean(trimmed), np.std(trimmed)
m, s

(np.float64(178.10278947124948), np.float64(7.017054887136004))

这是一个代表具有与修剪数据相同的均值和标准差的正态分布的 NormalCdf 对象。

from  thinkstats  import NormalCdf

cdf_normal_model = NormalCdf(m, s, name='normal model')

这是它与数据 CDF 的比较。

cdf_height = Cdf.from_seq(heights, name="data")
cdf_normal_model.plot(ls=":", color="gray")
cdf_height.step()

xlim = [140, 210]
decorate(xlabel="Height (cm)", ylabel="CDF", xlim=xlim)

使用 gaussian_kde 来创建一个近似男性身高的 PDF。提示：调查 bw_method 参数，它可以用来控制估计密度的平滑度。绘制估计密度，并将其与均值为 m 和标准差为 s 的 NormalPdf 进行比较。

使用 NormalPdf 和 area_under 来计算在正常模型中身高在 178 至 185 厘米之间的人的比例。使用 NormalCdf 来计算相同比例，并检查结果是否一致。最后，使用数据的经验 Cdf 来查看数据集中有多少人在同一范围内。

Think Stats: Exploratory Data Analysis in Python, 3rd Edition

代码许可：MIT 许可证

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 国际

变量之间的关系

原文：allendowney.github.io/ThinkStats/chap07.html

到目前为止，我们一次只查看一个变量。在本章中，我们开始研究变量之间的关系。如果知道一个变量就能给你关于另一个变量的信息，那么这两个变量就是相关的。例如，身高和体重是相关的——个子高的人往往比较重。当然，这并不是一个完美的关系：也有又矮又重的人和又高又轻的人。但如果你试图猜测某人的体重，知道他们的身高会比不知道更准确。

本章介绍了多种可视化变量之间关系的方法，以及一种量化关系强度，即相关性的方法。

点击此处运行此笔记本在 Colab 上。

from  os.path  import basename, exists

def  download(url):
    filename = basename(url)
    if not exists(filename):
        from  urllib.request  import urlretrieve

        local, _ = urlretrieve(url, filename)
        print("Downloaded " + local)

download("https://github.com/AllenDowney/ThinkStats/raw/v3/nb/thinkstats.py")

try:
    import  empiricaldist
except ImportError:
    %pip install empiricaldist

import  numpy  as  np
import  pandas  as  pd
import  matplotlib.pyplot  as  plt

from  thinkstats  import decorate

散点图

如果你遇到一个数学特别擅长的人，你会期望他们的语言技能比平均水平更好还是更差？一方面，你可能会想象人们会在一个或另一个领域专长，所以擅长一个领域的人可能在另一个领域不那么擅长。另一方面，你可能会期望一个一般都很聪明的人在两个领域都高于平均水平。让我们找出哪个是正确的。

我们将使用 1997 年国家青年纵向调查（NLSY97）的数据，该调查“追踪了 1980-84 年间出生的 8,984 名美国青年的生活”。公共数据集包括参与者在几个标准化测试中的得分，包括在大学入学中最常用的 SAT 和 ACT 测试。因为测试者会得到数学和语言部分的单独分数，我们可以使用这些数据来探索数学能力和语言能力之间的关系。

我使用 NLS Investigator 创建了一个摘录，其中包含我将用于此分析的自变量。在他们的许可下，我可以重新分发这个摘录。关于下载数据的说明在本书的笔记本中。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/nlsy97-extract.csv.gz")

我们可以使用read_csv来读取数据，并使用replace将缺失数据的特殊代码替换为np.nan。

missing_codes = [-1, -2, -3, -4, -5]
nlsy = pd.read_csv("nlsy97-extract.csv.gz").replace(missing_codes, np.nan)
nlsy.shape

(8984, 34)

nlsy.head()

| | R0000100 | R0490200 | R0536300 | R0536401 | R0536402 | R1235800 | R1318200 | R1482600 | R3961900 | R3989200 | ... | R9872200 | R9872300 | R9872400 | S1552700 | U0008900 | U1845500 | U3444000 | U4949700 | Z9083800 | Z9083900 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | 1 | NaN | 2 | 9 | 1981 | 1 | NaN | 4 | NaN | NaN | ... | 293.0 | 250.0 | 333.0 | NaN | 120000.0 | NaN | NaN | NaN | 16.0 | 4.0 |
| 1 | 2 | NaN | 1 | 7 | 1982 | 1 | 145.0 | 2 | NaN | NaN | ... | 114.0 | 230.0 | 143.0 | NaN | 98928.0 | 116000.0 | 188857.0 | 180000.0 | 14.0 | 2.0 |
| 2 | 3 | NaN | 2 | 9 | 1983 | 1 | 82.0 | 2 | NaN | NaN | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | 75000.0 | 16.0 | 4.0 |
| 3 | 4 | NaN | 2 | 2 | 1981 | 1 | NaN | 2 | NaN | NaN | ... | 195.0 | 230.0 | 216.0 | NaN | 85000.0 | 45000.0 | NaN | NaN | 13.0 | 2.0 |
| 4 | 5 | NaN | 1 | 10 | 1982 | 1 | NaN | 2 | NaN | NaN | ... | 293.0 | 230.0 | 231.0 | NaN | 210000.0 | 212000.0 | NaN | 240000.0 | 12.0 | 2.0 |

5 行 × 34 列

DataFrame包含调查中每个 8984 名参与者的一个行，以及我选择的每个 34 个变量的一个列。列名本身没有太多意义，所以让我们用我们将使用的更可解释的名称来替换它们。

nlsy["sat_verbal"] = nlsy["R9793800"]
nlsy["sat_math"] = nlsy["R9793900"]

两个列都包含一些小于 200 的值，这是不可能的，因为 200 是最低分，所以我们将它们替换为np.nan。

columns = ["sat_verbal", "sat_math"]

for column in columns:
    invalid = nlsy[column] < 200
    nlsy.loc[invalid, column] = np.nan

接下来，我们将使用dropna来选择只有两个分数都有效的行。

nlsy_valid = nlsy.dropna(subset=columns).copy()
nlsy_valid.shape

(1398, 36)

SAT 分数是标准化的，所以平均分是 500，标准差是 100。在 NLSY 样本中，平均值和标准差接近这些值。

sat_verbal = nlsy_valid["sat_verbal"]
sat_verbal.mean(), sat_verbal.std()

(np.float64(501.80972818311875), np.float64(108.36562024213643))

sat_math = nlsy_valid["sat_math"]
sat_math.mean(), sat_math.std()

(np.float64(503.0829756795422), np.float64(109.8329973731453))

现在，为了看看这些变量之间是否存在关系，让我们看看一个散点图。

plt.scatter(sat_verbal, sat_math)

decorate(xlabel="SAT Verbal", ylabel="SAT Math")

使用scatter函数的默认选项，我们可以看到关系的一般形状。在测试的一个部分表现好的人，在其他部分也倾向于表现更好。

然而，这个版本的图形是过度绘图的，这意味着有很多重叠的点，这可能会产生误导性的关系印象。点密度最高的中心区域并不像应该的那样暗——相比之下，极端值比应该的更暗。过度绘图往往会给异常值过多的视觉权重。

我们可以通过减小标记的大小来改进图形，这样它们就重叠得更少。

plt.scatter(sat_verbal, sat_math, s=5)

decorate(xlabel="SAT Verbal", ylabel="SAT Math")

现在我们可以看到标记已经按行和列对齐，因为分数被四舍五入到最接近的 10 的倍数。在这个过程中，一些信息丢失了。

我们无法恢复这些信息，但我们可以通过抖动数据来最小化对散点图的影响，这意味着向数据中添加随机噪声以抵消四舍五入的影响。以下函数接受一个序列，并通过添加具有均值为 0 和给定标准差的正态分布的随机值来抖动它。结果是 NumPy 数组。

def  jitter(seq, std=1):
    n = len(seq)
    return np.random.normal(0, std, n) + seq

如果我们用标准差为 3 的分数进行抖动，散点图中的行和列就不再可见了。

sat_verbal_jittered = jitter(sat_verbal, 3)
sat_math_jittered = jitter(sat_math, 3)

plt.scatter(sat_verbal_jittered, sat_math_jittered, s=5)

decorate(xlabel="SAT Verbal", ylabel="SAT Math")

抖动减少了四舍五入的视觉影响，使关系的形状更清晰。但一般来说，你应该只为可视化目的抖动数据，避免使用抖动数据进行分析。

在这个例子中，即使调整了标记大小并抖动数据，仍然存在一些重叠。所以让我们再试一次：我们可以使用alpha关键字使标记部分透明。

plt.scatter(sat_verbal_jittered, sat_math_jittered, s=5, alpha=0.2)

decorate(xlabel="SAT Verbal", ylabel="SAT Math")

在透明度下，重叠的数据点看起来更暗，因此暗度与密度成正比。

虽然散点图是一种简单且广泛使用的可视化方法，但要正确使用它们可能很困难。一般来说，需要一些尝试和错误来调整标记大小、透明度和抖动，以找到变量之间关系的最佳视觉表示。

十分位数图

散点图提供了变量之间关系的总体印象，但还有其他可视化方法可以更深入地了解关系的本质。其中之一是十分位数图。

要生成十分位数图，我们将受访者按语言分数排序，并将他们分成 10 组，称为十分位数。我们可以使用qcut方法计算十分位数。

deciles = pd.qcut(nlsy_valid["sat_verbal"], 10, labels=False) + 1
deciles.value_counts().sort_index()

sat_verbal
1     142
2     150
3     139
4     140
5     159
6     130
7     148
8     121
9     138
10    131
Name: count, dtype: int64

每个十分位数的受访者数量大致相等。

现在，我们可以使用groupby方法根据decile将DataFrame分成组。

df_groupby = nlsy_valid.groupby(deciles)
df_groupby

<pandas.core.groupby.generic.DataFrameGroupBy object at 0x7fa0adfd49d0>

结果是一个DataFrameGroupBy对象，表示组。我们可以从中选择sat_math列。

series_groupby = df_groupby["sat_math"]
series_groupby

<pandas.core.groupby.generic.SeriesGroupBy object at 0x7fa0b01af7d0>

结果是一个SeriesGroupBy对象，表示每个十分位数的数学分数。我们可以使用quantile函数计算每个组的第 10、50 和 90 个百分位数。

low = series_groupby.quantile(0.1)
median = series_groupby.quantile(0.5)
high = series_groupby.quantile(0.9)

十分位数图显示了每个十分位数组的这些百分位数。在下面的图中，线表示中位数，阴影区域表示第 10 和第 90 百分位数之间的区域。

xs = median.index
plt.fill_between(xs, low, high, alpha=0.2)
plt.plot(xs, median, label="median")

decorate(xlabel="SAT Verbal Decile", ylabel="SAT Math")

作为一种替代方案，我们可以计算每个小组的中位数语言分数，并将这些值绘制在 x 轴上，而不是十分位数数字。

xs = df_groupby["sat_verbal"].median()

plt.fill_between(xs, low, high, alpha=0.2)
plt.plot(xs, median, color="C0", label="median")

decorate(xlabel="SAT Verbal", ylabel="SAT Math")

看起来这些变量之间的关系是线性的——也就是说，中位数语言分数的每次增加都对应着中位数数学分数的大致相等增加。

更普遍地说，我们可以将受访者分成任意数量的组，不一定非要是 10 个，我们可以在每个组中计算其他汇总统计量，而不仅仅是这些百分位数。

排名相关性

NLSY 是纵向的，这意味着它随着时间的推移跟踪同一组人。我们一直在研究的群体包括 1980 年至 1984 年出生的人。那些参加 SAT 考试的人可能在 1990 年代末参加，当时他们大约 18 岁。因此，当他们在 2021 年谈论他们的收入时，他们已经 30 多岁或 40 多岁。让我们给收入数据列一个更易理解的名称。

nlsy["income"] = nlsy["U4949700"]
nlsy["income"].describe()

count      6051.000000
mean     104274.239960
std      108470.571497
min           0.000000
25%       38000.000000
50%       80000.000000
75%      134157.000000
max      599728.000000
Name: income, dtype: float64

这一列的值是毛家庭收入，即受访者及其家庭其他成员从所有来源报告的总收入，以美元（USD）为单位。以下是收入分布的图示。

cdf_income = Cdf.from_seq(nlsy["income"])
cdf_income.step()

decorate(xlabel="Income (USD)", ylabel="CDF")

注意到接近 600,000 的步骤——超过此阈值的值被限制以保护参与者的匿名性。现在，这是受访者 SAT 数学成绩和其未来收入的散点图。

scatter(nlsy, "piat_math", "income")

decorate(xlabel="PIAT math", ylabel="Gross Family Income (USD)")

看起来这些变量之间存在关系。以下是相关系数。

corrcoef(nlsy, "piat_math", "income")

np.float64(0.30338587288641233)

相关系数约为 0.3，这意味着如果某人在 15 岁时 PIAT 数学成绩比平均值高出一个标准差，我们预计他们在 40 岁时收入将比平均值高出约 0.3 个标准差。这不如 PIAT 成绩和 SAT 成绩之间的相关性强，但考虑到影响收入的因素数量，这相当强。

实际上，皮尔逊相关系数可能低估了关系的强度。正如我们在之前的散点图中可以看到的，两个变量在极端值方面都有明显的过剩。因为相关系数基于与平均值的偏差的乘积，它对这些极端值很敏感。

一个更稳健的替代方法是秩相关，它基于分数的秩而不是标准化分数。我们可以使用 Pandas 的rank方法来计算每个分数和每个收入的秩。

valid = nlsy.dropna(subset=["piat_math", "income"])

piat_math_rank = valid["piat_math"].rank(method="first")
income_rank = valid["income"].rank(method="first")

使用method="first"参数，rank将秩从 1 分配到序列的长度，即 4101。

income_rank.min(), income_rank.max()

(np.float64(1.0), np.float64(4101.0))

这是一张收入秩与数学分数秩的散点图。

plt.scatter(piat_math_rank, income_rank, s=5, alpha=0.2)

decorate(xlabel="PIAT math rank", ylabel="Income rank")

这是秩之间的相关性。

np.corrcoef(piat_math_rank, income_rank)[0, 1]

np.float64(0.38148396696764847)

结果大约是 0.38，略高于皮尔逊相关系数，后者为 0.30。因为秩相关对极端值的影响不太敏感，所以它可能是衡量这些变量之间关系强度的一个更好的指标。

thinkplot提供了一个封装了本节代码的rankcorr函数。

from  thinkstats  import rankcorr

rankcorr(nlsy, "piat_math", "income")

np.float64(0.38474681505344815)

SciPy 提供了一个类似的功能，称为spearmanr，因为秩相关也称为斯皮尔曼相关。

from  scipy.stats  import spearmanr

spearmanr(valid["piat_math"], valid["income"]).statistic

np.float64(0.38474681505344815)

作为练习，你将有机会使用皮尔逊相关和秩相关来计算 SAT 语文分数和收入之间的相关性。

术语表

散点图：一种可视化，通过为数据集中的每个观测值绘制一个点来显示两个变量之间的关系。
过度绘图：如果许多标记重叠，散点图就会过度绘图，这使得难以区分不同密度的区域，可能会错误地表示关系。
抖动：在图表中的数据点上添加随机噪声，以使重叠的值更易于可见。
十分位图：一种将数据根据一个变量分成十分位（十个组）的图表，然后对每个组总结另一个变量。
十分位：通过排序数据并将其分成大约相等的十部分所创建的组之一。
皮尔逊相关系数：一个统计量，用于衡量两个变量之间线性关系的强度和符号（正或负）。
标准分数：一个经过标准化处理的量，表示它相对于平均值的标准差。
相关矩阵：一个表格，显示数据集中每对变量的相关系数。
秩相关：一种通过使用值的秩而不是实际值来量化关系强度的稳健方法。
随机对照试验：一种实验，其中受试者被随机分配到接受不同治疗的组。
干预组：在实验中，接受正在测试的干预的组。
对照组：在实验中，不接受干预或接受已知效果的治疗的组。
自然实验：一种使用自然形成的组进行的实验，有时可以模仿随机分配。
因果推断：识别和量化因果关系的方法。

练习

练习 7.1

thinkstats模块提供了一个名为decile_plot的函数，它封装了本章前面的代码。我们可以这样调用它来可视化 SAT 语文和数学分数之间的关系。

from  thinkstats  import decile_plot

decile_plot(nlsy, "sat_verbal", "sat_math")
decorate(xlabel="SAT Verbal", ylabel="SAT Math")

制作 PIAT 数学分数和收入的十分位图。这似乎是线性关系吗？

练习 7.2

绘制收入与 SAT 数学分数的散点图。计算皮尔逊相关系数和等级相关系数。它们有显著差异吗？

绘制收入与 SAT 语文分数的散点图，并计算两个相关系数。哪个是未来收入的更强预测指标？

练习 7.3

让我们看看学生的高中平均成绩（GPA）与他们的 SAT 分数之间的相关性。这是 NLSY 数据集中编码 GPA 的变量。

missing_codes = [-6, -7, -8, -9]
nlsy["gpa"] = nlsy["R9871900"].replace(missing_codes, np.nan) / 100
nlsy["gpa"].describe()

count    6004.000000
mean        2.818408
std         0.616357
min         0.100000
25%         2.430000
50%         2.860000
75%         3.260000
max         4.170000
Name: gpa, dtype: float64

这就是 GPA 分布看起来像什么。

cdf_income = Cdf.from_seq(nlsy["gpa"])
cdf_income.step()
decorate(xlabel="GPA", ylabel="CDF")

绘制 GPA 与 SAT 数学分数之间的关系散点图，并计算相关系数。对于 GPA 与 SAT 语文分数之间的关系也做同样处理。哪个 SAT 分数是 GPA 更好的预测指标？

练习 7.4

让我们调查教育和收入之间的关系。NLSY 数据集包括一个列，报告了每个受访者的最高学位。这些值被编码为整数。

nlsy["degree"] = nlsy["Z9083900"]
nlsy["degree"].value_counts().sort_index()

degree
0.0     877
1.0    1167
2.0    3531
3.0     766
4.0    1713
5.0     704
6.0      64
7.0     130
Name: count, dtype: int64

但我们可以使用这些列表来解码它们。

positions = [0, 1, 2, 3, 4, 5, 6, 7]
labels = [
    "None",
    "GED",
    "High school diploma",
    "Associate's degree",
    "Bachelor's degree",
    "Master's degree",
    "PhD",
    "Professional degree",
]

并制作一个表示教育成就分布的Pmf。

from  empiricaldist  import Pmf

Pmf.from_seq(nlsy["degree"]).bar()

plt.xticks(positions, labels, rotation=30, ha="right")
decorate(ylabel="PMF")

绘制income与degree的散点图。为了避免重叠，对degree的值进行抖动，并调整标记大小和透明度。

使用groupby方法按degree对受访者进行分组。从DataFrameGroupBy对象中选择income列；然后使用quantile方法计算每个组的均值、10 分位数和 90 分位数。使用fill_between绘制 10 分位数和 90 分位数之间的区域，并使用plot绘制均值。

你能说些什么关于每增加一个学位所关联的收入溢价？

练习 7.4

行为风险因素监测系统（BRFSS）数据集包括大约 40 万受访者的自报身高和体重。下载数据的说明在本章的笔记本中。

绘制身高与体重的散点图。你可能需要抖动数据以模糊由于四舍五入而可见的行和列。由于样本量很大，你将不得不调整标记大小和透明度以避免重叠。此外，由于两种测量都有异常值，你可能想使用xlim和ylim来放大覆盖大多数受访者的区域。

这是我们如何加载数据的方法。

download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/CDBRFS08.ASC.gz")

from  thinkstats  import read_brfss

brfss = read_brfss()
brfss["htm3"].describe()

count    409129.000000
mean        168.825190
std          10.352653
min          61.000000
25%         160.000000
50%         168.000000
75%         175.000000
max         236.000000
Name: htm3, dtype: float64

制作体重与身高的十分位图。这种关系看起来是线性的吗？计算相关系数和等级相关系数。它们有显著差异吗？你认为哪一个更好地量化了这些变量之间的关系？

Think Stats: Exploratory Data Analysis in Python, 3rd Edition

代码许可：MIT License

文本许可：Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International

posted @ 2025-10-02 10:20 绝不原创的飞龙阅读(13) 评论(0) 收藏举报

刷新页面返回顶部

龙哥盟

掠夺·扩张·投机·博弈

tnkstat3e-merge-1

统计思维（程序员的概率统计）第三版（二）

模型分布

二项分布

泊松分布

指数分布

正态分布

为什么建模？

术语表

练习

练习 5.1

练习 5.2

练习 5.3

概率密度函数

比较分布

概率密度

指数概率密度函数（PDF）

比较概率质量函数（PMFs）和概率密度函数（PDFs）

核密度估计

分布框架

术语表

练习

练习 6.1

练习 6.2

变量之间的关系

散点图

十分位数图

相关系数

相关系数的强度

排名相关性

相关性与因果关系

术语表

练习

练习 7.1

练习 7.2

练习 7.3

练习 7.4

练习 7.4

公告