tnkstat3e-merge-3

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

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

译者：飞龙

协议：CC BY-NC-SA 4.0

时间序列分析

allendowney.github.io/ThinkStats/chap12.html

时间序列是从随时间变化的系统中获取的一系列测量值。我们在前几章中使用的大多数工具，如回归，也可以用于时间序列。但还有一些特别适用于此类数据的方法。

作为例子，我们将查看两个数据集：2001 年至 2024 年间美国的可再生能源发电量，以及同一时间段的天气数据。我们将开发将时间序列分解为长期趋势和重复季节性成分的方法。我们将使用线性回归模型来拟合和预测趋势。我们还将尝试一个广泛使用的分析时间序列数据模型，正式名称为“自回归积分移动平均”，简称 ARIMA。

Think Stats 的第三版现在可以从 Bookshop.org 和 Amazon 购买（这些是联盟链接）。如果你正在享受免费的在线版本，考虑 为我买杯咖啡。

点击此处运行此笔记本在 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

plt.rcParams["figure.dpi"] = 300 

电力

作为时间序列数据的例子，我们将使用美国能源信息署的数据集——它包括 2001 年至 2024 年间来自可再生能源的每月总发电量。下载数据的说明在本书的笔记本中。

下一个单元格将下载数据，我于 2024 年 9 月 17 日从 www.eia.gov/electricity/data/browser/ 下载了这些数据。

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

加载数据后，我们必须进行一些转换，以便将其转换为易于处理的形式。

elec = (
    pd.read_csv("Net_generation_for_all_sectors.csv", skiprows=4)
    .drop(columns=["units", "source key"])
    .set_index("description")
    .replace("--", np.nan)
    .transpose()
    .astype(float)
) 

在重新格式化的数据集中，每一列都是每月总发电量的序列，单位为千兆瓦时（GWh）。以下是列标签，显示了不同的电力来源或“部门”。

elec.columns 

Index(['Net generation for all sectors', 'United States',
       'United States : all fuels (utility-scale)', 'United States : nuclear',
       'United States : conventional hydroelectric',
       'United States : other renewables', 'United States : wind',
       'United States : all utility-scale solar', 'United States : geothermal',
       'United States : biomass',
       'United States : hydro-electric pumped storage',
       'United States : all solar',
       'United States : small-scale solar photovoltaic'],
      dtype='object', name='description') 

索引中的标签是表示月份和年份的字符串——以下是前 12 个月份。

elec.index[:12] 

Index(['Jan 2001', 'Feb 2001', 'Mar 2001', 'Apr 2001', 'May 2001', 'Jun 2001',
       'Jul 2001', 'Aug 2001', 'Sep 2001', 'Oct 2001', 'Nov 2001', 'Dec 2001'],
      dtype='object') 

如果我们将这些字符串替换为 Pandas Timestamp 对象，将更容易处理这些数据。我们可以使用 date_range 函数生成一系列 Timestamp 对象，从 2001 年 1 月开始，频率代码为 "ME"，代表“月末”，因此填充了每个月的最后一天。

elec.index = pd.date_range(start="2001-01", periods=len(elec), freq="ME")
elec.index[:6] 

DatetimeIndex(['2001-01-31', '2001-02-28', '2001-03-31', '2001-04-30',
               '2001-05-31', '2001-06-30'],
              dtype='datetime64[ns]', freq='ME') 

现在索引是一个 DataTimeIndex，数据类型为 datetime64[ns]，这是在 NumPy 中定义的——64 表示每个标签使用 64 位，ns 表示它具有纳秒精度。

分解

作为第一个例子，我们将看看从 2001 年 1 月到 2024 年 6 月的核电站发电量是如何变化的，并将时间序列分解为长期趋势和周期成分。以下是美国核电站的月度发电总量。

actual_options = dict(color="C0", lw=1, alpha=0.6)
trend_options = dict(color="C1", alpha=0.6)
pred_options = dict(color="C2", alpha=0.6, ls=':')
model_options = dict(color="gray", alpha=0.6, ls='--') 

nuclear = elec["United States : nuclear"]
nuclear.plot(label="nuclear", **actual_options)

decorate(ylabel="GWh") 

看起来有一些增加和减少，但它们很难清楚地看到，因为每个月之间有很大的变化。为了更清楚地看到长期趋势，我们可以使用 rolling 和 mean 方法来计算移动平均。

trend = nuclear.rolling(window=12).mean() 

window=12 参数选择重叠的 12 个月份的区间，因此第一个区间包含从第一个开始的 12 个测量值，第二个区间包含从第二个开始的 12 个测量值，依此类推。对于每个区间，我们计算平均产量。

下面是结果的样子，包括原始数据。

nuclear.plot(label="nuclear", **actual_options)
trend.plot(label="trend", **trend_options)
decorate(ylabel="GWh") 

趋势仍然相当多变。我们可以通过使用更长的窗口来进一步平滑它，但现在我们将坚持使用 12 个月的窗口。

如果我们从原始数据中减去趋势，结果是一个“去趋势”的时间序列，这意味着长期平均值接近常数。下面是它的样子。

detrended = (nuclear - trend).dropna()
detrended.plot(label="detrended", **actual_options)
decorate(ylabel="GWh") 

看起来存在一个重复的年度模式，这是有道理的，因为电力的需求从季节到季节是变化的，它被用来在冬天产生热量和在夏天运行空调。为了描述这个年度模式，我们可以选择索引中 datetime 对象的月份部分，按月份对数据进行分组，并计算平均产量。下面是月平均产量的样子。

monthly_averages = detrended.groupby(detrended.index.month).mean()
monthly_averages.plot(label="monthly average", **actual_options)
decorate(ylabel="GWh") 

在 x 轴上，第 1 个月是 1 月，第 12 个月是 12 月。电力生产在最冷和最热的月份最高，在 4 月和 10 月最低。

我们可以使用 monthly_averages 来构建数据的季节成分，这是一个与 nuclear 长度相同的序列，其中每个月的元素是该月的平均值。下面是它的样子。

seasonal = monthly_averages[nuclear.index.month]
seasonal.index = nuclear.index
seasonal.plot(label="seasonal", **actual_options)
decorate(ylabel="GWh") 

每个 12 个月份的周期与其他周期相同。

趋势成分和季节成分的总和代表每个月的预期值。

expected = trend + seasonal 

下面是它与原始序列相比的样子。

expected.plot(label="expected", **pred_options)
nuclear.plot(label="actual", **actual_options)
decorate(ylabel="GWh") 

如果我们从原始序列中减去这个总和，结果就是残差成分，它代表了每个月相对于预期值的偏离。

resid = nuclear - expected
resid.plot(label="residual", **actual_options)
decorate(ylabel="GWh") 

我们可以将残差视为影响能源生产的所有因素的总和，但这些因素既不是长期趋势也不是季节性成分所解释的。其中包括天气、因维护而停机的设备，以及由于特定事件引起的需求变化。由于残差是许多不可预测的，有时甚至无法知晓的因素的总和，我们通常将其视为一个随机量。

这就是残差的分布情况。

from  thinkstats  import plot_kde

plot_kde(resid.dropna())
decorate(xlabel="Residual (GWh)", ylabel="Density") 

它类似于正态分布的钟形曲线，这与它是许多随机贡献的总和的假设是一致的。

为了量化这个模型如何描述原始序列，我们可以计算确定系数，它表示残差方差相对于原始序列方差的大小。

rsquared = 1 - resid.var() / nuclear.var()
rsquared 

np.float64(0.9054559977517084) 

(R²)值约为 0.92，这意味着长期趋势和季节性成分解释了序列中 92%的变化。这个(R²)值比我们之前章节中看到的要高得多，但在时间序列数据中这是常见的——尤其是在我们构建的模型类似于数据的这种情况下。

我们刚才所经历的过程被称为季节性分解。StatsModels 提供了一个执行此操作的函数，称为seasonal_decompose。

from  statsmodels.tsa.seasonal  import seasonal_decompose

decomposition = seasonal_decompose(nuclear, model="additive", period=12) 

model="additive"参数表示加法模型，因此序列被分解为趋势、季节性成分和残差的总和。我们很快就会看到乘法模型。period=12参数表示季节性成分的持续时间为 12 个月。

结果是一个包含三个成分的对象。本章的笔记本提供了一个绘制它们的函数。

def  plot_decomposition(original, decomposition):
    plt.figure(figsize=(6, 5))

    ax1 = plt.subplot(4, 1, 1)
    plt.plot(original, label="Original", color="C0", lw=1)
    plt.ylabel("Original")

    plt.subplot(4, 1, 2, sharex=ax1)
    plt.plot(decomposition.trend, label="Trend", color="C1", lw=1)
    plt.ylabel("Trend")

    plt.subplot(4, 1, 3, sharex=ax1)
    plt.plot(decomposition.seasonal, label="Seasonal", color="C2", lw=1)
    plt.ylabel("Seasonal")

    plt.subplot(4, 1, 4, sharex=ax1)
    plt.plot(decomposition.resid, label="Residual", color="C3", lw=1)
    plt.ylabel("Residual")

    plt.tight_layout() 

plot_decomposition(nuclear, decomposition) 

结果与我们自己计算的结果相似，但由于实现细节的不同，存在一些小的差异。

这种季节性分解有助于了解时间序列的结构。正如我们将在下一节中看到的，它对于进行预测也是很有用的。

预测

我们可以使用季节分解的结果来预测未来。为了演示，我们将使用以下函数将时间序列拆分为一个训练序列，我们将用它来生成预测，以及一个测试序列，我们将用它来查看它们是否准确。

def  split_series(series, n=60):
    training = series.iloc[:-n]
    test = series.iloc[-n:]
    return training, test 

在 n=60 的情况下，测试序列的持续时间为五年，始于 2019 年 7 月。

training, test = split_series(nuclear)
test.index[0] 

Timestamp('2019-07-31 00:00:00') 

现在，假设是 2019 年 6 月，你被要求生成核发电厂五年电力生产的预测。为了回答这个问题，我们将使用训练数据来构建模型，然后使用该模型来生成预测。我们将从训练数据的季节分解开始。

decomposition = seasonal_decompose(training, model="additive", period=12)
trend = decomposition.trend 

现在，我们将对趋势拟合一个线性模型。解释变量 months 是从序列开始到当前月份的月份数。

import  statsmodels.formula.api  as  smf

months = np.arange(len(trend))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months", data=data).fit() 

这是结果总结。

from  thinkstats  import display_summary

display_summary(results) 

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	6.482e+04	131.524	492.869	0.000	6.46e+04	6.51e+04
months	10.9886	1.044	10.530	0.000	8.931	13.046
R-squared:	0.3477

(R²) 值约为 0.35，这表明模型并不特别适合数据。我们可以通过绘制拟合线来更好地理解这一点。我们将使用 predict 方法来计算训练数据和测试数据的预期值。

months = np.arange(len(training) + len(test))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = nuclear.index 

这是趋势成分和线性模型。

trend.plot(**trend_options)
pred_trend.plot(label="linear model", **model_options)
decorate(ylabel="GWh") 

有很多事情没有被线性模型捕捉到，但看起来总体上有一个上升趋势。

接下来，我们将使用分解中的季节成分来计算一系列月平均数。

seasonal = decomposition.seasonal
monthly_averages = seasonal.groupby(seasonal.index.month).mean() 

我们可以通过在 monthly_averages 中查找拟合线的日期来预测季节成分。

pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index 

最后，为了生成预测，我们将季节成分添加到趋势中。

pred = pred_trend + pred_seasonal 

这是训练数据和预测值。

pred.plot(label="prediction", **pred_options)
training.plot(label="training", **actual_options)
decorate(ylabel="GWh") 

预测值与训练数据拟合得相当好，基于长期趋势将继续的假设，预测看起来像是一个合理的预测。

现在，从未来的角度来看，让我们看看这个预测有多准确。以下是 2019 年 7 月至五年间的预测值和实际值。

forecast = pred[test.index]
forecast.plot(label="predicted", **pred_options)
test.plot(label="actual", **actual_options)
decorate(ylabel="GWh") 

预测的第一年相当不错，但 2020 年核反应堆的生产低于预期——可能是由于 COVID-19 大流行——并且它从未回到长期趋势。

为了量化预测的准确性，我们将使用平均绝对百分比误差（MAPE），该函数如下计算。

def  MAPE(predicted, actual):
    ape = np.abs(predicted - actual) / actual
    return np.mean(ape) * 100 

在这个例子中，预测的平均误差为 3.81%。

MAPE(forecast, test) 

np.float64(3.811940747879257) 

我们将在本章的后面回到这个例子，看看我们是否可以用不同的模型做得更好。

乘法模型

在上一节中我们使用的加法模型假设时间序列是长期趋势、季节性成分和残差的和，这意味着季节性成分和残差的幅度不会随时间变化。

作为违反这个假设的一个例子，让我们看看自 2014 年以来小规模太阳能发电的情况。

solar = elec["United States : small-scale solar photovoltaic"].dropna()
solar.plot(label="solar", **actual_options)
decorate(ylabel="GWh") 

在这个区间内，总产量增加了几倍。很明显，季节性变化的幅度也增加了。

如果我们假设季节性和随机变化的幅度与趋势的幅度成正比，那么这表明了一种替代加法模型的方法，其中时间序列是三个成分的乘积。

要尝试这个乘法模型，我们将这个序列分为训练集和测试集。

training, test = split_series(solar) 

使用seasonal_decompose函数并传入model="multiplicative"参数。

decomposition = seasonal_decompose(training, model="multiplicative", period=12) 

这就是结果看起来像什么。

plot_decomposition(training, decomposition) 

现在的季节性和残差成分是乘法因子。因此，季节性成分似乎从低于趋势约 25%变化到高于趋势约 25%。残差成分通常小于 5%，除了第一个周期中的一些较大因子。我们可以这样提取模型成分。

trend = decomposition.trend
seasonal = decomposition.seasonal
resid = decomposition.resid 

这个模型的(R²)值非常高。

rsquared = 1 - resid.var() / training.var()
rsquared 

np.float64(0.9999999992978134) 

太阳能板的产量很大程度上取决于它所暴露的阳光，因此产量遵循年度周期是有意义的。

为了预测长期趋势，我们将使用二次模型。

months = range(len(training))
data = pd.DataFrame({"trend": trend, "months": months}).dropna()
results = smf.ols("trend ~ months + I(months**2)", data=data).fit() 

在 Patsy 公式中，子串I(months**2)向模型添加了一个二次项，所以我们不必显式地计算它。以下是结果。

display_summary(results) 

	系数	标准误差	t	P>|t|	[0.025	0.975]
截距	766.1962	13.494	56.782	0.000	739.106	793.286
月份	22.2153	0.938	23.673	0.000	20.331	24.099
I(months ** 2)	0.1762	0.014	12.480	0.000	0.148	0.205
R-squared:	0.9983

线性和二次项的 p 值非常小，这表明二次模型比线性模型能捕捉到更多关于趋势的信息——并且(R²)值非常高。

现在我们可以使用这个模型来计算过去和未来的趋势的预期值。

months = range(len(solar))
df = pd.DataFrame({"months": months})
pred_trend = results.predict(df)
pred_trend.index = solar.index 

这就是它的样子。

pred_trend.plot(label="quadratic model", **model_options)
trend.plot(**trend_options)
decorate(ylabel="GWh") 

二次模型很好地拟合了过去趋势。现在我们可以使用季节性成分来预测未来的季节性变化。

monthly_averages = seasonal.groupby(seasonal.index.month).mean()
pred_seasonal = monthly_averages[pred_trend.index.month]
pred_seasonal.index = pred_trend.index 

最后，为了计算过去值的回溯预测和未来的预测，我们将趋势和季节成分相乘。

pred = pred_trend * pred_seasonal 

这里是结果以及训练数据。

training.plot(label="training", **actual_options)
pred.plot(label="prediction", **pred_options)
decorate(ylabel="GWh") 

这些预测与训练数据拟合得很好，预测看起来是合理的——现在让我们看看它们是否最终是准确的。以下是预测和测试数据。

future = pred[test.index]
future.plot(label="prediction", **pred_options)
test.plot(label="actual", **actual_options)
decorate(ylabel="GWh") 

在前三年，预测非常好。之后，看起来实际增长超过了预期。

在这个例子中，季节分解对于建模和预测太阳能生产效果很好，但在先前的例子中，它对核能生产并不非常有效。在下一节中，我们将尝试不同的方法，即自回归。

自回归

自回归的第一个想法是未来将与过去相似。例如，在我们迄今为止查看的时间序列中，有一个明显的年度周期。所以如果你被要求预测下个月的预测，一个好的起点是上个月。

为了看看这可能有多好，让我们回到nuclear，它包含核发电厂的月度电力生产，并计算连续年份中同一月份的差异，这些差异被称为“年对年”差异。

diff = (nuclear - nuclear.shift(12)).dropna()
diff.plot(label="year over year differences", **actual_options)
decorate(ylabel="GWh") 

这些差异的幅度远小于原始序列的幅度，这表明自回归的第二个想法，即预测这些差异可能比预测原始值更容易。

为了达到这个目的，让我们看看差异序列中连续元素之间是否存在相关性。如果有的话，我们可以使用这些相关性根据先前值预测未来值。

我将从创建一个DataFrame开始，将差异放在第一列，并将相同差异（滞后 1 个月、2 个月和 3 个月）放入后续列。这些列被命名为lag1、lag2和lag3，因为它们包含的序列已经被滞后或延迟。

df_ar = pd.DataFrame({"diff": diff})
for lag in [1, 2, 3]:
    df_ar[f"lag{lag}"] = diff.shift(lag)

df_ar = df_ar.dropna() 

这里是这些列之间的相关性。

df_ar.corr()[["diff"]] 

	diff
diff	1.000000
lag1	0.562212
lag2	0.292454
lag3	0.222228

这些相关性被称为滞后相关性或自相关——前缀“auto”表示我们正在取序列与其自身的相关性。作为一个特殊情况，diff和lag1之间的相关性被称为序列相关性，因为它是在序列中连续元素之间的相关性。

这些相关性足够强，足以表明它们应该有助于预测，所以让我们将它们放入多元回归中。以下函数使用DataFrame的列来创建一个 Patsy 公式，其中第一列作为响应变量，其他列作为解释变量。

def  make_formula(df):
  """Make a Patsy formula from column names."""
    y = df.columns[0]
    xs = " + ".join(df.columns[1:])
    return f"{y} ~ {xs}" 

这里是一个线性模型的结果，该模型根据前三个值预测序列中的下一个值。

formula = make_formula(df_ar)
results_ar = smf.ols(formula=formula, data=df_ar).fit()
display_summary(results_ar) 

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	24.2674	114.674	0.212	0.833	-201.528	250.063
lag1	0.5847	0.061	9.528	0.000	0.464	0.706
lag2	-0.0908	0.071	-1.277	0.203	-0.231	0.049
lag3	0.1026	0.062	1.666	0.097	-0.019	0.224
R-squared:	0.3239

现在，我们可以使用predict方法生成序列中过去值的预测。以下是这些回溯预测与数据相比的图示。

pred_ar = results_ar.predict(df_ar)
pred_ar.plot(label="predictions", **pred_options)
diff.plot(label="differences", **actual_options)
decorate(ylabel="GWh") 

预测在某些地方很好，但(R²)值只有大约 0.319，所以还有改进的空间。

resid_ar = (diff - pred_ar).dropna()
R2 = 1 - resid_ar.var() / diff.var()
R2 

np.float64(0.3190252265690783) 

提高预测的一种方法是从该模型计算残差，并使用另一个模型来预测残差——这是自回归的第三个想法。

移动平均

假设现在是 2019 年 6 月，有人要求你预测 2020 年 6 月的值。你的第一个猜测可能是今年的值将在明年重复。

现在假设是 2020 年 5 月，有人要求你修改你对 2020 年 6 月的预测。你可以使用过去三个月的结果和上一节中的自相关模型来预测年同比增长率。

最后，假设你检查了最近几个月的预测，并发现它们一直过低。这表明下个月的预测也可能过低，因此你可以将其向上调整。基本假设是最近的预测误差可以预测未来的预测误差。

为了看看它们是否真的有效，我们可以创建一个DataFrame，其中第一列包含自回归模型的残差，其他列包含残差的滞后版本。在这个例子中，我将使用 1 个月和 6 个月的滞后。

df_ma = pd.DataFrame({"resid": resid_ar})

for lag in [1, 6]:
    df_ma[f"lag{lag}"] = resid_ar.shift(lag)

df_ma = df_ma.dropna() 

我们可以使用ols为残差创建一个自回归模型。这个模型部分被称为“移动平均”，因为它以类似于移动平均效果的方式减少了预测中的变异性。我发现这个术语并不特别有帮助，但它是一种惯例。

无论如何，这里是对残差自回归模型的总结。

formula = make_formula(df_ma)
results_ma = smf.ols(formula=formula, data=df_ma).fit()
display_summary(results_ma) 

	coef	std err	t	P>|t|	[0.025	0.975]
Intercept	-14.0016	114.697	-0.122	0.903	-239.863	211.860
lag1	0.0014	0.062	0.023	0.982	-0.120	0.123
lag6	-0.1592	0.063	-2.547	0.011	-0.282	-0.036
R-squared:	0.0247

(R²) 值相当小，所以看起来这个模型的部分不会很有帮助。但 6 个月滞后值的 p 值很小，这表明它比我们预期的偶然性贡献了更多信息。

现在，我们可以使用该模型来生成残差的反预测。

pred_ma = results_ma.predict(df_ma) 

然后，为了生成逐年差异的反预测，我们将第二个模型的调整值添加到第一个模型的反预测中。

pred_diff = pred_ar + pred_ma 

两个模型的总和的 (R²) 值约为 0.332，这略好于没有移动平均调整的结果（0.319）。

resid_ma = (diff - pred_diff).dropna()
R2 = 1 - resid_ma.var() / diff.var()
R2 

np.float64(0.3315101001391231) 

接下来，我们将利用这些逐年差异来生成原始值的反预测。

基于自回归的反预测

要生成反预测，我们首先将逐年差异放入与原始索引对齐的 Series 中。

pred_diff = pd.Series(pred_diff, index=nuclear.index) 

使用 isna 检查 NaN 值，我们发现新 Series 的前 21 个元素缺失。

n_missing = pred_diff.isna().sum()
n_missing 

np.int64(21) 

这是因为我们通过 12 个月将 Series 平移来计算逐年差异，然后我们将第一个自回归模型的差异平移了 3 个月，并将第一个模型的残差平移了 6 个月用于第二个模型。每次我们像这样平移一个 Series，我们都会在开头丢失一些值，这些平移的总和是 21。

因此，在我们能够生成反预测之前，我们必须通过将原始数据的前 21 个元素复制到一个新的 Series 中来预热泵。

pred_series = pd.Series(index=nuclear.index, dtype=float)
pred_series.iloc[:n_missing] = nuclear.iloc[:n_missing] 

现在，我们可以运行以下循环，它从索引 21（即第 22 个元素）填充到末尾。每个元素是前一年的值和预测的逐年差异的总和。

for i in range(n_missing, len(pred_series)):
    pred_series.iloc[i] = pred_series.iloc[i - 12] + pred_diff.iloc[i] 

现在，我们将用 NaN 替换我们复制的元素，这样我们就不因“完美预测”前 21 个值而获得信用。

pred_series[:n_missing] = np.nan 

这是与原始值相比的反预测的样子。

pred_series.plot(label="predicted", **pred_options)
nuclear.plot(label="actual", **actual_options)
decorate(ylabel="GWh") 

它们看起来相当不错，(R²) 值约为 0.86。

resid = (nuclear - pred_series).dropna()
R2 = 1 - resid.var() / nuclear.var()
R2 

np.float64(0.8586566911201015) 

我们用来计算这些反预测的模型称为 SARIMA，它是被称为 ARIMA 的一组模型之一。这些缩写的每个部分都指代模型的一个元素。

• S 代表季节性，因为第一步是计算相隔一个季节周期的值之间的差异。

• AR 代表自回归，我们用它来建模差异中的滞后相关性。

• I 代表积分，因为我们用来计算 pred_series 的迭代过程与微积分中的积分类似。

• MA 代表移动平均，这是我们在第一个模型的残差上运行的第二个自回归模型的常规名称。

ARIMA 模型是建模时间序列数据的强大且多功能的工具。

ARIMA

StatsModel 提供了一个名为 tsa 的库，代表“时间序列分析”——它包括一个名为 ARIMA 的函数，用于拟合 ARIMA 模型并生成预测。

为了拟合我们在前几节中开发的 SARIMA 模型，我们将使用两个元组作为参数调用此函数：order 和 seasonal_order。以下是 order 中与我们在前几节中使用的模型相对应的值。

order = ([1, 2, 3], 0, [1, 6]) 

order 中的值表示：

• 应该在 AR 模型中包含哪些滞后项——在这个例子中是前三个。

• 它应该计算连续元素之间差异的次数——在这个例子中是 0，因为我们计算了季节差分，我们将在下一分钟讨论这一点。

• 应该在 MA 模型中包含哪些滞后项——在这个例子中是第一个和第六个。

现在是 seasonal_order 中的值。

seasonal_order = (0, 1, 0, 12) 

第一个和第三个元素是 0，这意味着该模型不包含季节性 AR 或季节性 MA。第二个元素是 1，这意味着它计算季节差分——最后一个元素是季节周期。

这是我们将如何使用 ARIMA 来制作和拟合这个模型的示例。

import  statsmodels.tsa.api  as  tsa

model = tsa.ARIMA(nuclear, order=order, seasonal_order=seasonal_order)
results_arima = model.fit()
display_summary(results_arima) 

	coef	std err	z	P>|z|	[0.025	0.975]
ar.L1	0.0458	0.379	0.121	0.904	-0.697	0.788
ar.L2	-0.0035	0.116	-0.030	0.976	-0.230	0.223
ar.L3	0.0375	0.049	0.769	0.442	-0.058	0.133
ma.L1	0.2154	0.382	0.564	0.573	-0.533	0.964
ma.L6	-0.0672	0.019	-3.500	0.000	-0.105	-0.030
sigma2	3.473e+06	1.9e-07	1.83e+13	0.000	3.47e+06	3.47e+06

结果包括 AR 模型中三个滞后项、MA 模型中两个滞后项以及 sigma2 的估计系数，sigma2 是残差的方差。

从 results_arima 中我们可以提取 fittedvalues，它包含回溯预测。由于我们在计算回溯预测时在开始处有缺失值，因此在 fittedvalues 的开始处也有不正确的值，我们将丢弃这些值。

fittedvalues = results_arima.fittedvalues[n_missing:] 

拟合值与我们计算的类似，但并不完全相同——可能是因为 ARIMA 处理初始条件的方式不同。

fittedvalues.plot(label="ARIMA model", **pred_options)
nuclear.plot(label="actual", **actual_options)
decorate(ylabel="GWh") 

(R²) 值也类似，但并不完全相同。

resid = fittedvalues - nuclear
R2 = 1 - resid.var() / nuclear.var()
R2 

np.float64(0.8262717330822065) 

ARIMA 函数使得实验不同的模型版本变得容易。

作为练习，尝试在 order 和 seasonal_order 中使用不同的值，看看你是否能找到一个具有更高 (R²) 的模型。

使用 ARIMA 进行预测

ARIMA 返回的对象提供了一个名为 get_forecast 的方法，该方法生成预测。为了演示，我们将时间序列分为训练集和测试集，并将相同的模型拟合到训练集。

training, test = split_series(nuclear)
model = tsa.ARIMA(training, order=order, seasonal_order=seasonal_order)
results_training = model.fit() 

我们可以使用结果来为测试集生成预测。

forecast = results_training.get_forecast(steps=len(test)) 

结果是一个包含一个名为 forecast_mean 的属性和一个返回置信区间的函数的对象。

forecast_mean = forecast.predicted_mean
forecast_ci = forecast.conf_int()
forecast_ci.columns = ["lower", "upper"] 

我们可以像这样绘制结果，并与实际的时间序列进行比较。

plt.fill_between(
    forecast_ci.index,
    forecast_ci.lower,
    forecast_ci.upper,
    lw=0,
    color="gray",
    alpha=0.2,
)
plt.plot(forecast_mean.index, forecast_mean, label="forecast", **pred_options)
plt.plot(test.index, test, label="actual", **actual_options)
decorate(ylabel="GWh") 

实际值几乎完全位于预测的置信区间内。以下是预测的 MAPE。

MAPE(forecast_mean, test) 

np.float64(3.3817549248044547) 

预测的平均误差为 3.38%，略好于我们从季节分解中得到的结果（3.81%）。

ARIMA 比季节分解更灵活，并且通常能做出更好的预测。在这个时间序列中，自相关性并不特别强，因此 ARIMA 的优势是适度的。

术语表

• 时间序列：一个数据集，其中每个值都与特定的时间相关联，通常表示在固定间隔内进行的测量。

• 季节分解：一种将时间序列分解为长期趋势、重复的季节成分和残差成分的方法。

• 训练序列：用于拟合模型的时间序列的一部分。

• 测试序列：用于检查模型生成的预测准确性的时间序列的一部分。

• 回溯预测：对过去观察到的值的预测，通常用于测试或验证模型。

• 窗口：时间序列中连续的值序列，用于计算移动平均。

• 移动平均：通过平均重叠窗口中的值来平滑波动的时间序列。

• 序列相关性：时间序列连续元素之间的相关性。

• 自相关性：时间序列与其自身位移或滞后版本之间的相关性。

• 滞后：序列相关性或自相关性中位移的大小。

练习

练习 12.1

作为季节分解的一个例子，让我们来模拟美国月平均地表温度。我们将使用来自 Our World in Data 的数据集，该数据集包括“从 1950 年到 2024 年，大多数国家地面以上 2 米处测量的空气温度[摄氏度]，包括陆地、海洋和内陆水面”，用于大多数国家。数据下载说明在本书的笔记本中。

# The following cell downloads data prepared by Our World in Data,
# which I Downloaded September 18, 2024
# from https://ourworldindata.org/grapher/average-monthly-surface-temperature

# Based on modified data from Copernicus Climate Change Service information (2019)
# with "major processing" by Our World in Data 

filename = "monthly-average-surface-temperatures-by-year.csv"
download("https://github.com/AllenDowney/ThinkStats/raw/v3/data/" + filename) 

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

temp = pd.read_csv("monthly-average-surface-temperatures-by-year.csv") 

temp.head() 

| | 实体 | 代码 | 年份 | 2024 | 2023 | 2022 | 2021 | 2020 | 2019 | 2018 | ... | 1959 | 1958 | 1956 | 1954 | 1952 | 1957 | 1955 | 1953 | 1951 | 1950 |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| 0 | 阿富汗 | AFG | 1 | 3.300064 | -4.335608 | -0.322859 | -1.001608 | -2.560545 | 0.585145 | 1.042471 | ... | -2.333814 | 0.576404 | -3.351925 | -2.276692 | -2.812619 | -4.239172 | -2.191683 | -2.915993 | -3.126317 | -2.655707 |
| 1 | 阿富汗 | AFG | 2 | 1.024550 | 4.187041 | 2.165870 | 5.688000 | 2.880046 | 0.068664 | 3.622793 | ... | -1.545529 | 0.264962 | 0.455350 | -0.304205 | 0.798226 | -2.747945 | 1.999074 | 1.983414 | -2.642800 | -3.996040 |
| 2 | 阿富汗 | AFG | 3 | 5.843506 | 10.105444 | 10.483686 | 9.777976 | 6.916731 | 5.758049 | 10.794412 | ... | 5.942937 | 7.716459 | 5.090270 | 4.357703 | 4.796146 | 4.434027 | 7.066073 | 4.590406 | 3.054388 | 3.491112 |
| 3 | 阿富汗 | AFG | 4 | 11.627398 | 14.277164 | 17.227650 | 15.168276 | 12.686832 | 13.838840 | 14.321226 | ... | 13.752827 | 14.712909 | 11.982360 | 12.155265 | 13.119270 | 8.263829 | 10.418768 | 11.087193 | 9.682878 | 8.332797 |
| 4 | 阿富汗 | AFG | 5 | 18.957850 | 19.078170 | 19.962734 | 19.885902 | 18.884047 | 18.461287 | 18.100782 | ... | 17.388723 | 16.352045 | 20.125462 | 18.432117 | 17.614851 | 15.505956 | 15.599709 | 17.865084 | 17.095737 | 17.329062 |

5 行 × 78 列

下面的单元格从 2001 年到系列结束选择美国的数据，并将其打包到一个 Pandas Series中。

temp_us = temp.query("Code == 'USA'")
columns = [str(year) for year in range(2000, 2025)]
temp_series = temp_us.loc[:, columns].transpose().stack()
temp_series.index = pd.date_range(start="2000-01", periods=len(temp_series), freq="ME") 

这就是它的样子。

temp_series.plot(label="monthly average", **actual_options)
decorate(ylabel="Surface temperature (degC)") 

毫不奇怪，这里有一个明显的季节性模式。使用 12 个月为周期的加性季节分解。将线性模型拟合到趋势线。在这个区间内，地表温度的平均年增长是多少？如果你好奇，可以用其他区间或来自其他国家的数据重复这个分析。

练习 12.2

在本章的早期，我们使用乘性季节分解来模拟 2014 年至 2019 年小规模太阳能发电量，并预测 2019 年至 2024 年的产量。现在让我们用同样的方法来处理规模化的太阳能发电。这就是时间序列的样子。

util_solar = elec["United States : all utility-scale solar"].dropna()
util_solar = util_solar[util_solar.index.year >= 2014]
util_solar.plot(**actual_options)
decorate(ylabel="GWh") 

使用split_series将数据分割成训练和测试系列。对训练系列进行 12 个月周期的乘性分解。将线性或二次模型拟合到趋势，并生成包括季节成分在内的五年预测。将预测与测试系列一起绘制，并计算平均绝对百分比误差（MAPE）。

练习 12.3

让我们看看 ARIMA 模型如何拟合美国水力发电站的产量。以下是 2001 年至 2024 年的时间序列图。

hydro = elec["United States : conventional hydroelectric"]
hydro.plot(**actual_options)
decorate(ylabel="GWh") 

使用 12 个月季节周期的 SARIMA 模型拟合这些数据。在模型的自回归和移动平均部分尝试不同的滞后，看看是否能找到一个最大化模型(R²)值的组合。生成五年预测并绘制其置信区间。

注意：根据你包含在模型中的滞后，你可能会发现拟合值的前 12 到 24 个元素不可靠。你可能想在绘图或计算(R²)之前移除它们。

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

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可证

文本许可：Creative Commons 知识共享署名-非商业性使用-相同方式共享 4.0 国际

生存分析

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

生存分析是一种描述事物持续时间的分析方法。它通常用于研究人类寿命，但也适用于机械和电子组件的“生存”，或者更普遍地，用于任何事件发生之前的时间间隔——甚至空间间隔。

我们将从简单的例子开始，即灯泡的寿命，然后考虑一个更复杂的例子，即首次结婚的年龄以及过去 50 年美国这一变化。

点击此处运行此笔记本在 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 

生存函数

生存分析中的一个基本概念是生存函数，它表示在给定持续时间后存活下来的群体比例。作为一个例子，我们将计算灯泡寿命的生存函数。

我们将使用 2007 年进行的一项实验的数据。研究人员安装了 50 个新的灯泡，并持续点亮。他们每 12 小时检查一次灯泡，并记录任何失效的灯泡的寿命——直到所有 50 个灯泡都失效。下载数据的说明在本书的笔记本中。

下面的单元格下载了数据，其文档在此链接中。

数据集来源：

V.J. Menon 和 D.C. Agrawal，灯丝灯泡的更新率：理论与实验。失效分析预防杂志。2007 年 12 月，第 421 页，表 2/ DOI：10.1007/s11668-007-9074-9

download(
    "https://gist.github.com/epogrebnyak/7933e16c0ad215742c4c104be4fbdeb1/raw/c932bc5b6aa6317770c4cbf43eb591511fec08f9/lamps.csv"
) 

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

df = pd.read_csv("lamps.csv", index_col=0)
df.tail() 

	h	f	K
i
---	---	---	---
28	1812	1	4
29	1836	1	3
30	1860	1	2
31	1980	1	1
32	2568	1	0

h列包含以小时为单位的寿命。f列记录了在h的每个值下失效的灯泡数量。为了表示寿命的分布，我们将这些值放入Pmf对象中并对其进行归一化。

from  empiricaldist  import Pmf

pmf_bulblife = Pmf(df["f"].values, index=df["h"])
pmf_bulblife.normalize() 

np.int64(50) 

我们可以使用make_cdf来计算累积分布函数（CDF），它表示在h的每个值或之前失效的灯泡的比例。例如，78%的灯泡在 1656 小时或之前失效。

cdf_bulblife = pmf_bulblife.make_cdf()
cdf_bulblife[1656] 

np.float64(0.7800000000000002) 

生存函数是指在h的每个值之后失效的灯泡的比例，它是累积分布函数的补数。因此，我们可以这样计算它。

complementary_cdf = 1 - cdf_bulblife
complementary_cdf[1656] 

np.float64(0.21999999999999975) 

22%的灯泡在 1656 小时后失效。

empiricaldist库提供了一个表示生存函数的Surv对象，以及一个名为make_surv的方法来创建一个生存函数。

surv_bulblife = cdf_bulblife.make_surv()
surv_bulblife[1656] 

np.float64(0.21999999999999997) 

如果我们绘制累积分布函数（CDF）和生存函数，我们可以看到它们是互补的——也就是说，在h的所有值上它们的和为 1。

cdf_bulblife.plot(ls="--", label="CDF")
surv_bulblife.plot(label="Survival")

decorate(xlabel="Light bulb duration (hours)", ylabel="Probability") 

在这个意义上，累积分布函数（CDF）和生存函数是等价的——如果我们给出其中一个，我们可以计算出另一个——但在生存分析中，更常见的是使用生存曲线。计算生存曲线是下一个重要概念，即失效函数的步骤。

失效函数

在灯泡数据集中，h的每个值代表一个以小时h结束的 12 小时间隔——我将称之为“间隔h”。假设我们知道一个灯泡已经存活到间隔h，我们想知道它在间隔h内失效的概率。为了回答这个问题，我们可以使用生存函数，它表示超过间隔h的灯泡的比例，以及 PMF，它表示在间隔h内失效的比例。这两个值的和是可能在间隔h内失效的灯泡的比例，这些灯泡被称为“处于风险中”。例如，26%的灯泡在 1656 间隔内处于风险中。

at_risk = pmf_bulblife + surv_bulblife
at_risk[1656] 

np.float64(0.25999999999999995) 

并且 4%的所有灯泡在 1656 间隔内失效。

pmf_bulblife[1656] 

np.float64(0.04) 

失效函数是pmf_bulblife和at_risk的比率。

hazard = pmf_bulblife / at_risk
hazard[1656] 

np.float64(0.15384615384615388) 

在所有存活到 1656 区间的灯泡中，大约有 15%在 1656 区间内失效。

我们不必自己计算失效函数，可以使用empiricaldist，它提供了一个表示失效函数的Hazard对象，以及一个计算它的make_hazard方法。

hazard_bulblife = surv_bulblife.make_hazard()
hazard_bulblife[1656] 

np.float64(0.15384615384615397) 

这里展示了灯泡的失效函数是什么样的。

hazard_bulblife.plot()
decorate(xlabel="Light bulb duration (hours)", ylabel="Hazard") 

我们可以看到在某些地方失效函数比其他地方高，但这种方式可视化失效函数可能会误导人，尤其是在我们数据不多的一部分范围内。一个更好的选择是绘制累积失效函数，它是失效函数的累积和。

cumulative_hazard = hazard_bulblife.cumsum()
cumulative_hazard.plot()

decorate(xlabel="Light bulb duration (hours)", ylabel="Cumulative hazard") 

在失效概率高的地方，累积失效函数是陡峭的。在失效概率低的地方，累积失效函数是平缓的。在这个例子中，我们可以看到在 1500 到 2000 小时之间失效函数最高。之后，失效函数下降——尽管这个结果仅基于一个寿命异常长的灯泡，所以在其他数据集中可能看起来不同。

既然我们已经了解了生存和失效函数的一般概念，让我们将它们应用到更大的数据集上。

婚姻数据

在许多国家，人们的结婚年龄比以前晚，而且越来越多的人保持未婚。为了探索美国这些趋势，我们将使用生存分析工具和国家家庭成长调查（NSFG）的数据。

我们在前面章节中使用的 NSFG 数据集是怀孕文件，它包含每个调查受访者报告的怀孕的行。在本章中，我们将处理受访者文件，其中包含有关受访者本人的信息。

我已经整理了 1982 年至 2019 年之间进行的九次调查的回应，并选择了与婚姻相关的数据。下载此摘录的说明在本书的笔记本中。

以下单元格下载了数据，这是一个 CSV 文件，我创建的，它结合了 1982 年至 2019 年 NSFG 调查的多个迭代的数据。

数据准备详情请参阅这个笔记本。

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

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

resp = pd.read_csv("marriage_nsfg_female.csv.gz")
resp.shape 

(70183, 34) 

摘录中包括超过 70,000 名受访者的每一行，并包含以下与年龄和婚姻相关的变量。

• cmbirth：受访者的出生日期，对所有受访者都已知。

• cmintvw：受访者被采访的日期，对所有受访者都已知。

• cmmarrhx：如果适用且已知，受访者的首次结婚日期。

• evrmarry：如果受访者在采访日期之前已经结婚，则为 1，否则为 0。

前三个变量编码为“世纪-月份”——即自 1899 年 12 月以来的整数月数。因此，世纪-月份 1 是 1900 年 1 月。

为了探索代际变化，我们将受访者按出生十年分组。我们将使用以下函数，它接受cmbirth的值并计算相应的出生十年。它使用整数除法运算符//除以 10 并向下取整。

month0 = pd.to_datetime("1899-12-31")

def  decade_of_birth(cmbirth):
    date = month0 + pd.DateOffset(months=cmbirth)
    return date.year // 10 * 10 

我们可以使用此函数和apply方法来计算每个受访者的出生十年，并将其分配给一个名为cohort的新列。在这个上下文中，群体是一群有共同点的人——比如他们出生的十年——在分析中作为一个群体被对待。

value_counts的结果显示了每个群体中的人数。

from  thinkstats  import value_counts

resp["cohort"] = resp["cmbirth"].apply(decade_of_birth)
value_counts(resp["cohort"]) 

cohort
1930      325
1940     3608
1950    10631
1960    14953
1970    16438
1980    14271
1990     8552
2000     1405
Name: count, dtype: int64 

数据集包括 1950 年代至 1980 年代每个十年出生的 10,000 多名受访者，以及早期和后期十年中受访者较少。

接下来，我们将计算每个受访者在结婚时的年龄（如果适用）以及他们在采访时的年龄。

resp["agemarr"] = (resp["cmmarrhx"] - resp["cmbirth"]) / 12
resp["age"] = (resp["cmintvw"] - resp["cmbirth"]) / 12 

要开始使用这些数据，我们将使用以下函数，它接受一个DataFrame和一个群体列表作为参数，并返回一个字典，将每个群体映射到一个Surv对象。对于每个群体，它选择他们首次结婚时的年龄，并使用Surv.from_seq来计算生存函数。dropna=False参数将NaN值包含在生存函数中，因此结果包括尚未结婚的人。

from  empiricaldist  import Surv

def  make_survival_map(resp, cohorts):
    surv_map = {}

    grouped = resp.groupby("cohort")
    for cohort in cohorts:
        group = grouped.get_group(cohort)
        surv_map[cohort] = Surv.from_seq(group["agemarr"], dropna=False)

    return surv_map 

这是使用此函数的方法。

cohorts = [1980, 1960, 1940]
surv_map = make_survival_map(resp, cohorts) 

下面是出生于 1940 年代、1960 年代和 1980 年代的人的结果。

import  matplotlib.pyplot  as  plt
import  cycler

# Extract the default color cycle
default_colors = plt.rcParams["axes.prop_cycle"].by_key()["color"]

# Define the desired line styles
linestyles = ["--", "-", "-.", "--", "-", "-.", "--", "-", "-.", "--"]

# Ensure we cycle through styles and colors properly
custom_cycler = cycler.cycler(color=default_colors) + cycler.cycler(
    linestyle=linestyles
)

# Apply the new cycler
plt.rc("axes", prop_cycle=custom_cycler) 

for cohort, surv in surv_map.items():
    surv.plot(label=f"{cohort}s")

ylim = [-0.05, 1.05]
decorate(xlabel="Age", ylabel="Prob never married", ylim=ylim) 

如果我们直接接受这些结果，它们表明早期世代的人结婚较早，而且更多的人最终结婚了。然而，我们不应立即解释这些结果，因为它们是不正确的。我们必须解决两个问题：

• 如第一章所述，NSFG 使用分层抽样，这意味着它故意对某些群体进行过抽样。

• 此外，这种计算生存函数的方法并没有充分考虑尚未结婚的人。

对于第一个问题，我们将使用一种称为加权自助法的重抽样方法。对于第二个问题，我们将使用一种称为 Kaplan-Meier 估计的方法。我们将从重抽样开始。

加权自助法

NSFG 数据集包括一个名为finalwgt的列，其中包含每个受访者的抽样权重，即他们所代表的总人口中的数量。我们可以在重抽样过程中使用这些权重来校正分层抽样。以下函数接受一个DataFrame和包含抽样权重的列名。它根据抽样权重重抽样DataFrame的行，并返回一个新的DataFrame。

def  resample_rows_weighted(df, column="finalwgt"):
    n = len(df)
    weights = df[column]
    return df.sample(n, weights=weights, replace=True) 

当前数据集包括来自调查几个迭代周期的受访者，称为周期，因此为了重抽样，我们必须按周期对受访者进行分组，对每个组进行重抽样，然后将组重新组合。这正是以下函数所做的事情。

def  resample_cycles(resp):
    grouped = resp.groupby("cycle")
    samples = [resample_rows_weighted(group) for _, group in grouped]
    return pd.concat(samples) 

要开始，我们将重抽样数据一次。

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

sample = resample_cycles(resp) 

之后我们将多次重抽样数据，以便我们可以看到由于随机抽样引起的多少变异。

下图显示了重抽样后的结果，与之前章节未重抽样的结果（以虚线表示）进行了比较。

for label, surv in surv_map.items():
    surv.plot(ls=":", color="gray", alpha=0.6)

survs_resampled = make_survival_map(sample, cohorts)

for label, surv in survs_resampled.items():
    surv.plot(label=label)

decorate(xlabel="Age", ylabel="Prob never married", ylim=ylim) 

plt.rc("axes", prop_cycle=plt.rcParamsDefault["axes.prop_cycle"]) 

有无重抽样的差异是显著的，这表明我们需要对分层抽样进行校正以获得准确的结果。

现在，让我们来解决第二个问题，处理不完整的数据。

估计风险函数

在灯泡示例中，我们知道所有 50 个灯泡的寿命，因此我们可以直接计算生存函数——我们还可以使用生存函数来计算风险函数。

在婚姻示例中，我们知道一些受访者在被采访之前已经结婚的初次结婚年龄。但对于从未结婚的受访者，我们不知道他们将来会在多大年龄结婚——或者是否会结婚。

这种缺失数据被称为截尾。这个术语可能听起来有些奇怪，因为截尾信息通常是故意隐藏的，但在这个情况下，它只是因为我们不知道未来。

然而，我们有一些可以工作的部分信息：如果某人在调查时未婚，我们知道他们结婚的年龄（如果他们结婚的话）必须大于他们当前的年龄。

我们可以使用这些部分信息来估计风险函数；然后我们可以使用风险函数来计算生存函数。这个过程被称为Kaplan-Meier 估计。

为了演示，我将从重新抽样的数据中选取一个队列。

resp60 = sample.query("cohort == 1960") 

对于在调查时已婚的受访者，我们将选择他们的初婚年龄。共有 9921 人，我们将它们称为“完整”案例。

complete = resp60.query("evrmarry == 1")["agemarr"]
complete.count() 

np.int64(9921) 

对于那些未婚的受访者，我们将选择他们在调查时的年龄。共有 5468 人，我们将它们称为“持续”案例。

ongoing = resp60.query("evrmarry == 0")["age"]
ongoing.count() 

np.int64(5468) 

现在，为了估计风险函数，我们将计算每个年龄的“处于风险中”的案例总数，包括到那个年龄为止所有未婚的人。创建一个FreqTab对象，该对象按年龄统计完整和持续案例的数量将很方便。

from  empiricaldist  import FreqTab

ft_complete = FreqTab.from_seq(complete)
ft_ongoing = FreqTab.from_seq(ongoing) 

例如，有 58 位受访者报告称他们在 25 岁时第一次结婚。

ft_complete[25] 

np.int64(58) 

还有 5 位在 25 岁时接受调查并报告称他们从未结婚的人。

ft_ongoing[25] 

np.int64(5) 

从这些FreqTab对象中，我们可以计算未归一化的Surv对象，这些对象包含超过每个年龄的完整和持续案例数量。

surv_complete = ft_complete.make_surv()
surv_ongoing = ft_ongoing.make_surv() 

例如，有 2848 人报告称他们在 25 岁之后结婚。

surv_complete[25] 

np.int64(2848) 

以及 25 岁之后接受调查且从未结婚的 2273 人。

surv_ongoing[25] 

np.int64(2273) 

我们刚刚计算的四个数字之和是处于风险中的受访者人数——也就是说，在 25 岁时可能结婚的人。术语“处于风险中”是医学中生存分析的遗产，它通常指疾病或死亡的风险。在婚姻的背景下，这似乎有些不一致，因为婚姻通常被认为是一个积极的里程碑。话虽如此，这是我们的计算方法。

at_risk = ft_complete[25] + ft_ongoing[25] + surv_complete[25] + surv_ongoing[25]
at_risk 

np.int64(5184) 

其中，实际在 25 岁结婚的人数是ft_complete[25]。因此，我们可以这样计算 25 岁时的风险函数。

hazard = ft_complete[25] / at_risk
hazard 

np.float64(0.011188271604938271) 

这就是我们如何计算单个年龄的风险函数。现在让我们计算整个函数，适用于所有年龄。我们将使用Index类的union方法来计算一个包含ft_complete和ft_ongoing中所有年龄的 Pandas Index。

ts = pd.Index.union(ft_complete.index, ft_ongoing.index) 

现在，我们可以通过在每个FreqTab和Surv对象中查找ts中的年龄来计算每个年龄处于风险中的人数。

at_risk = ft_complete(ts) + ft_ongoing(ts) + surv_complete(ts) + surv_ongoing(ts) 

最后，我们可以计算每个年龄的风险函数，并将结果放入Hazard对象中。

from  empiricaldist  import Hazard

hs = ft_complete(ts) / at_risk
hazard = Hazard(hs, ts) 

这就是累积风险函数的形状。

hazard.cumsum().plot()

decorate(xlabel="Age", ylabel="Cumulative hazard") 

它在 20 岁到 30 岁之间最陡峭，这意味着未婚者在这些年龄段的结婚“风险”最大。之后，累积风险水平稳定，这意味着风险逐渐降低。

估计生存函数

如果我们给定一个生存函数，我们知道如何计算风险函数。现在让我们反过来。

这里有一个思考方式。风险函数表示如果你尚未结婚，在每一个年龄结婚的概率。所以风险函数的补数是每个年龄保持未婚的概率。

为了“生存”到给定的年龄t，你必须在每个年龄保持未婚，包括t。做到这一点的概率是互补风险函数的乘积，我们可以这样计算。

ps = (1 - hazard).cumprod() 

Hazard对象有一个make_surv方法来完成这个计算。

surv = hazard.make_surv() 

下面是结果的样子，与之前的结果（虚线）相比，之前的结果纠正了分层重采样，但没有处理截尾数据。

survs_resampled[1960].plot(ls=":", color="gray", label="resampled")
surv.plot(label="Kaplan-Meier")

decorate(xlabel="Age", ylabel="Prob never married", ylim=ylim) 

我们可以看到正确处理截尾数据的重要性。

这样的生存函数是 1986 年一篇著名杂志文章的基础 – 《新闻周刊》报道称，40 岁的未婚女性“被恐怖分子杀害的可能性”比结婚的可能性更大。这个说法被广泛报道，并成为流行文化的一部分，但当时是错误的（因为它是基于错误的分析），后来证明更加错误（因为文化变化已经在进行中）。2006 年，《新闻周刊》又发表了一篇文章，承认他们的错误。

我鼓励你阅读更多关于这篇文章、其依据的统计数据以及公众的反应。这应该会提醒你，进行统计分析时要有道德上的责任，用适当的怀疑态度来解释结果，并准确、诚实地向公众展示。

Snopes 有关于这个主题的好文章。

以下函数封装了 Kaplan-Meier 估计的步骤。它接受完整和持续病例的生存时间序列作为参数，并返回一个Hazard对象。

def  estimate_hazard(complete, ongoing):
  """Kaplan-Meier estimation."""
    ft_complete = FreqTab.from_seq(complete)
    ft_ongoing = FreqTab.from_seq(ongoing)

    surv_complete = ft_complete.make_surv()
    surv_ongoing = ft_ongoing.make_surv()

    ts = pd.Index.union(ft_complete.index, ft_ongoing.index)
    at_risk = ft_complete(ts) + ft_ongoing(ts) + surv_complete(ts) + surv_ongoing(ts)

    hs = ft_complete(ts) / at_risk
    return Hazard(hs, ts) 

下面是一个函数，它接受一组受访者，提取生存时间，调用estimate_hazard来获取风险函数，然后计算相应的生存函数。

def  estimate_survival(group):
  """Estimate the survival function."""
    complete = group.query("evrmarry == 1")["agemarr"]
    ongoing = group.query("evrmarry == 0")["age"]
    hf = estimate_hazard(complete, ongoing)
    sf = hf.make_surv()
    return sf 

很快，我们将使用这些函数来计算生存函数的置信区间。但首先，让我们看看另一种计算 Kaplan-Meier 估计的方法。

Lifelines

一个名为lifelines的 Python 包提供了生存分析的工具，包括计算 Kaplan-Meier 估计的函数。

以下单元格在需要时安装lifelines。

try:
    import  lifelines
except ImportError:
    %pip install lifelines 

我们可以使用它来确认上一节中的结果是正确的。首先，我们将使用estimate_survival计算生存函数。

surv = estimate_survival(resp60) 

接下来，我们将使用lifelines来计算它。首先，我们需要将数据转换成lifelines所需的格式。

complete = complete.dropna()
durations = np.concatenate([complete, ongoing])
event_observed = np.concatenate([np.ones(len(complete)), np.zeros(len(ongoing))]) 

现在，我们可以创建一个KaplanMeierFitter对象并拟合数据。

from  lifelines  import KaplanMeierFitter

kmf = KaplanMeierFitter()
kmf.fit(durations=durations, event_observed=event_observed) 

<lifelines.KaplanMeierFitter:"KM_estimate", fitted with 15389 total observations, 5468 right-censored observations> 

在拟合数据后，我们可以调用plot函数来显示结果，包括估计的生存函数和置信区间——尽管在这个案例中置信区间是不正确的，因为它没有对分层抽样进行校正。

kmf.plot()

decorate(xlabel="Age", ylabel="Prob never married", ylim=ylim) 

与我们计算的生存函数不同，lifelines中的生存函数从 0 开始。但其余的函数在浮点误差范围内是相同的。

ps = kmf.survival_function_["KM_estimate"].drop(0)
np.allclose(ps, surv) 

True 

在下一节中，我们将使用加权重采样来计算考虑分层抽样的置信区间。

置信区间

我们计算的 Kaplan-Meier 估计是基于数据集的单次重采样。为了了解由于随机抽样引起的变异程度，我们将使用多次重采样进行分析，并绘制结果。

我们将使用以下函数，该函数接受一个DataFrame和一个队列列表，为每个队列估计生存函数，并返回一个字典，将每个整数队列映射到一个Surv对象。

此函数与make_survival_map相同，除了它调用estimate_survival，该函数使用 Kaplan-Meier 估计，而不是Surv.from_seq，后者仅在无删失数据的情况下工作。

def  estimate_survival_map(resp, cohorts):
  """Make a dictionary from cohorts to Surv objects."""
    surv_map = {}

    grouped = resp.groupby("cohort")
    for cohort in cohorts:
        group = grouped.get_group(cohort)
        surv_map[cohort] = estimate_survival(group)

    return surv_map 

以下循环生成了数据集的 101 次随机重采样，并创建了一个包含 101 个字典的列表，这些字典包含了估计的生存函数。

cohorts = [1940, 1950, 1960, 1970, 1980, 1990]

surv_maps = [estimate_survival_map(resample_cycles(resp), cohorts) for i in range(101)] 

为了绘制结果，我们将使用以下函数，该函数接受一个字典列表、一个整数队列和一个颜色字符串。它遍历字典，选择给定队列的生存函数，并以几乎透明的线条绘制它——这是可视化重采样之间变异的一种方法。

def  plot_cohort(surv_maps, cohort, color):
  """Plot results for a single cohort."""
    survs = [surv_map[cohort] for surv_map in surv_maps]
    for surv in survs:
        surv.plot(color=color, alpha=0.05)

    x, y = surv.index[-1], surv.iloc[-1]
    plt.text(x + 1, y, f"{cohort}s", ha="left", va="center") 

这里是 1940 年代到 1990 年代出生队列的结果。

colors = [f"C{i}" for i in range(len(cohorts))]

for cohort, color in zip(cohorts, colors):
    plot_cohort(surv_maps, cohort, color)

xlim = [8, 55]
decorate(xlabel="Age", ylabel="Prob never married", xlim=xlim, ylim=ylim) 

这种可视化对于探索来说已经足够好了，但线条看起来有些模糊，一些标签也重叠了。可能需要更多的工作来制作一个可用于发表的图表，但现在我们会保持简单。

可以看到几个明显的模式：

• 1940 年代出生的女性结婚最早——1950 年代和 1960 年代出生的队列结婚较晚，但未婚的比例大致相同。

• 1970 年代出生的女性结婚较晚，并且未婚的比例高于之前的队列。

• 1980 年代和 1990 年代出生的队列结婚的时间更晚，并且未婚的比例更高——尽管这些模式在未来可能会改变。

我们将不得不等待 NSFG 的下一个数据发布来了解更多信息。

预期剩余寿命

给定一个分布，我们可以计算预期剩余寿命作为已过时间的函数。例如，给定怀孕长度的分布，我们可以计算分娩的预期时间。为了演示，我们将使用 NSFG 的怀孕数据。

以下单元格下载数据文件并安装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 

我们将使用get_nsfg_groups来读取数据，并将其分为第一胎和其他。

from  nsfg  import get_nsfg_groups

live, firsts, others = get_nsfg_groups() 

我们将首先对数据进行一次重采样。

sample = resample_rows_weighted(live, "finalwgt") 

这是怀孕持续时间的 PMF。

pmf_durations = Pmf.from_seq(sample["prglngth"]) 

现在假设是怀孕的第 36 周的开始。记住最常见的怀孕长度是 39 周，我们预计剩余时间为 3-4 周。为了使这个估计更精确，我们可以识别分布中等于或超过 36 周的那些值。

t = 36
is_remaining = pmf_durations.qs >= t 

接下来，我们将创建一个新的Pmf对象，它只包含那些值，并将它们左移，使当前时间为 0。

ps = pmf_durations.ps[is_remaining]
qs = pmf_durations.qs[is_remaining] - t

pmf_remaining = Pmf(ps, qs) 

由于我们选择了Pmf中的值的一个子集，概率不再加起来为 1，但我们可以归一化Pmf，使它们加起来为 1。

pmf_remaining.normalize() 

np.float64(0.9155006558810669) 

这是结果，它显示了第 36 周开始时的剩余时间分布。

pmf_remaining.bar(label="Week 36")
decorate(xlabel="Remaining duration (weeks)", ylabel="PMF") 

这个分布的均值是预期的剩余时间。

pmf_remaining.mean() 

np.float64(3.2145671641791043) 

以下函数封装了这些步骤，并计算给定时间t下Pmf的剩余时间分布。

def  compute_pmf_remaining(pmf, t):
  """Distribution of remaining time."""
    is_remaining = pmf.qs >= t
    ps = pmf.ps[is_remaining]
    qs = pmf.qs[is_remaining] - t
    pmf_remaining = Pmf(ps, qs)
    pmf_remaining.normalize()
    return pmf_remaining 

以下函数接受一个怀孕长度的Pmf，并计算从第 36 周到第 43 周每周开始的预期剩余时间。

def  expected_remaining(pmf):
    index = range(36, 44)
    expected = pd.Series(index=index)

    for t in index:
        pmf_remaining = compute_pmf_remaining(pmf, t)
        expected[t] = pmf_remaining.mean()

    return expected 

这是数据单次重采样的结果。

expected = expected_remaining(pmf_durations)
expected 

36    3.214567
37    2.337714
38    1.479095
39    0.610133
40    0.912517
41    0.784211
42    0.582301
43    0.589372
dtype: float64 

为了看到由于随机采样引起的多少变化，我们可以运行多次重采样分析，并绘制结果。

for i in range(21):
    sample = resample_rows_weighted(live, "finalwgt")
    pmf_durations = Pmf.from_seq(sample["prglngth"])
    expected = expected_remaining(pmf_durations)
    expected.plot(color="C0", alpha=0.1)

decorate(
    xlabel="Weeks of pregnancy", ylabel="Expected remaining time (weeks)", ylim=[0, 3.4]
) 

在第 36 周到第 39 周之间，预期的剩余时间逐渐减少，直到在第 39 周开始时，大约是 0.6 周。但之后，曲线趋于平稳。在第 40 周开始时，预期的剩余时间仍然接近 0.6 周——实际上略高——在第 41 周、第 42 周和第 43 周开始时，几乎相同。对于焦急等待婴儿出生的人来说，这种行为似乎相当残酷。

术语表

• 生存分析：一组描述和预测感兴趣事件发生时间的的方法，通常关注寿命或持续时间。

• 生存函数：一个将时间t映射到超过t的生存概率的函数。

• 风险函数：一个将t映射到在t时刻经历事件的案例比例的函数，这些案例在t时刻之前幸存。

• 累积风险函数：风险函数的累积和，通常对可视化很有用。

• 加权 bootstrap：一种使用采样权重来纠正分层抽样的重采样形式，通过模拟代表性样本。

• 被审查的数据：仅部分已知的数据，因为感兴趣的事件尚未发生或未被观察到。

• Kaplan-Meier 估计：一种用于估计具有截尾观察值的数据集中的生存和风险函数的方法。

• 队列：具有共享特征（如出生十年）的一组主题，作为一个组进行分析。

练习

练习 13.1

我们可以使用本章中的方法估计婚姻持续时间的风险和生存函数。为了使事情简单，我们将只考虑第一次婚姻，并且我们将关注离婚作为终点，而不是分离或死亡。

在 NSFG 数据中，cmdivorcx 列包含每个受访者的第一次婚姻的离婚日期（如果适用），以世纪-月份编码。计算以离婚结束的婚姻的持续时间，以及至今仍在进行的婚姻的持续时间。

• 对于完整案例，计算 cmdivorcx 和 cmmarrhx 之间的经过时间。如果这两个值都有效——不是 NaN ——这意味着受访者的第一次婚姻以离婚结束。

• 为了识别持续案例，选择只结过一次婚并且仍然已婚的人。你可以使用 fmarno，它记录了每个受访者结婚的次数，以及 fmarital，它编码了他们的婚姻状况——值 1 表示受访者已婚。

在某些情况下，这些变量的值只是近似值，因此你可能会发现一小部分负数差异，但它们不应超过一年。

估计婚姻持续时间的风险和生存函数。绘制累积风险函数——离婚的危险性最高是在什么时候？绘制生存函数——有多少比例的婚姻以离婚结束？

# I suggest you use a single resampling of the data

sample = resample_cycles(resp) 

练习 13.2

2012 年，南加州大学的一组人口统计学家估计了 19 世纪初和 20 世纪初出生在瑞典的人的预期寿命。对于 0 至 91 岁的年龄，他们估计了年龄特定的死亡率，即给定年龄死亡的人数占所有在该年龄之前存活下来的人的比例——你可能认识这是风险函数。

我使用在线图形数字化工具从他们的论文中的图表获取数据，并将其存储在 CSV 文件中。下载数据的说明在本书的笔记本中。

数据来源：Beltrán-Sánchez, H., Crimmins, E. M., & Finch, C. E. (2012). 早期队列死亡率预测队列的衰老速度：历史分析。发育起源的健康与疾病杂志，3(5)，380-386。

下面的单元格下载了数据。

download(
    "https://github.com/AllenDowney/ThinkStats/raw/v3/data/mortality_rates_beltran2012.csv"
) 

我们可以像这样加载数据。

mortality = pd.read_csv("mortality_rates_beltran2012.csv", header=[0, 1]).dropna() 

下面的函数将数据插值以创建一个具有从 0 到 99 岁每个年龄的近似死亡率的危险函数。

from  scipy.interpolate  import interp1d
from  empiricaldist  import Hazard

def  make_hazard(ages, rates):
  """Make a Hazard function by interpolating a Series.

 series: Series

 returns: Hazard
 """
    interp = interp1d(ages, rates, fill_value="extrapolate")
    xs = np.arange(0, 100)
    ys = np.exp(interp(xs))
    return Hazard(ys, xs) 

现在我们可以创建一个这样的 Hazard 对象。

ages = mortality["1800", "X"].values
rates = mortality["1800", "Y"].values
hazard = make_hazard(ages, rates) 

这就是死亡率的样子。

hazard.plot()

decorate(xlabel="Age (years)", ylabel="Hazard") 

使用 make_surv 来根据这些比率制作生存函数，并使用 make_cdf 来计算相应的累积分布函数（CDF）。绘制结果。

然后使用 make_pmf 来制作表示寿命分布的 Pmf 对象，并绘制它。最后，使用 compute_pmf_remaining 来计算从 0 岁到 99 岁每个年龄的平均剩余寿命。绘制结果。

在剩余寿命曲线上，你应该会看到一个反直觉的模式——在生命的最初几年，剩余寿命会增加。因为 19 世纪初期的婴儿死亡率非常高，所以预期较大的孩子比年幼的孩子活得更久。大约 5 岁之后，预期寿命会回到我们预期的模式——年轻人预期比老年人活得更久。

如果你对这个主题感兴趣，可能会喜欢我的书《很可能是在过度思考》的第五章，该章节展示了来自统计学许多领域的类似反直觉结果。

可从 probablyoverthinking.it 获取。

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

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可证

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

解析方法

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

本书一直关注计算方法，如模拟和重采样，但我们解决的一些问题有解析解，这些解可以更快地计算。

本章介绍了一些这些方法，并解释了它们是如何工作的。在章节末尾，我提出了一些建议，关于如何将计算和分析方法整合到数据分析中。

在此处运行此笔记本。

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 

正态概率图

许多解析方法基于正态分布的性质，原因有两个：现实世界中许多测量的分布很好地近似于正态分布，而且正态分布具有使它们对分析有用的数学性质。

为了证明第一个观点，我们将查看企鹅数据集中的某些测量值。然后我们将探索正态分布的数学特性。下载数据的说明在本书的笔记本中。

下一个单元格将从 Allison Horst 创建的存储库中下载数据。

Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: 帕默群岛（南极）企鹅数据。R 包版本 0.1.0。allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218.

数据是作为导致这篇论文的研究的一部分收集的：Gorman KB, Williams TD, Fraser WR (2014).南极企鹅（属 Pygoscelis）群落中的生态性两性差异和环境可变性。PLoS ONE 9(3):e90081。doi.org/10.1371/journal.pone.0090081

download(
    "https://raw.githubusercontent.com/allisonhorst/palmerpenguins/c19a904462482430170bfe2c718775ddb7dbb885/inst/extdata/penguins_raw.csv"
) 

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

penguins = pd.read_csv("penguins_raw.csv")
penguins.shape 

(344, 17) 

数据集包含三种企鹅物种的测量值。在这个例子中，我们将选择阿德利企鹅。

adelie = penguins.query('Species.str.startswith("Adelie")').copy()
len(adelie) 

152 

为了查看企鹅体重是否遵循正态分布，我们将计算数据的经验累积分布函数（CDF）。

from  empiricaldist  import Cdf

weights = adelie["Body Mass (g)"].dropna()
cdf_weights = Cdf.from_seq(weights) 

我们将计算具有相同均值和标准差的正态分布的解析累积分布函数（CDF）。

m, s = weights.mean(), weights.std()
m, s 

(np.float64(3700.662251655629), np.float64(458.5661259101348)) 

from  scipy.stats  import norm

dist = norm(m, s)
qs = np.linspace(m - 3.5 * s, m + 3.5 * s)
ps = dist.cdf(qs) 

这是数据与正态模型相比的累积分布函数（CDF）的样子。

model_options = dict(color="gray", alpha=0.5, label="model")
plt.plot(qs, ps, **model_options)
cdf_weights.plot(label="data")

decorate(ylabel="CDF") 

正态分布可能足以作为这些数据的模型，但它肯定不是完美的拟合。

通常，绘制数据的累积分布函数（CDF）和模型的累积分布函数是一种评估模型拟合数据好坏的好方法。但这种方法的一个缺点是它依赖于我们对模型参数估计的好坏——在这个例子中，是均值和标准差。

另一个选择是正态概率图，它不依赖于我们估计参数的能力。在正态概率图中，(y)值是排序后的测量值。

ys = np.sort(weights) 

而(x)值是正态分布的相应百分位数，使用norm对象的ppf方法计算，该方法计算“百分位数函数”，即逆 CDF。

n = len(weights)
ps = (np.arange(n) + 0.5) / n
xs = norm.ppf(ps) 

如果测量值实际上是从正态分布中抽取的，那么(y)和(x)值应该落在一条直线上。为了看看它们做得如何，我们可以使用linregress来拟合一条线。

from  scipy.stats  import linregress

results = linregress(xs, ys)
intercept, slope = results.intercept, results.slope

fit_xs = np.linspace(-3, 3)
fit_ys = intercept + slope * fit_xs 

下图显示了(x)和(y)值以及拟合的线。

plt.plot(fit_xs, fit_ys, **model_options)
plt.plot(xs, ys, label="data")

decorate(xlabel="Standard normal", ylabel="Body mass (g)") 

正态概率图不是一条完美的直线，这表明正态分布不是这个数据的完美模型。

原因之一是数据集包括雄性和雌性企鹅，两组的平均值不同——让我们看看如果我们分别绘制这两组会发生什么。以下函数封装了我们用来制作正态概率图的步骤。

def  normal_probability_plot(sample, **options):
  """Makes a normal probability plot with a fitted line."""
    n = len(sample)
    ps = (np.arange(n) + 0.5) / n
    xs = norm.ppf(ps)
    ys = np.sort(sample)

    results = linregress(xs, ys)
    intercept, slope = results.intercept, results.slope

    fit_xs = np.linspace(-3, 3)
    fit_ys = intercept + slope * fit_xs

    plt.plot(fit_xs, fit_ys, color="gray", alpha=0.5)
    plt.plot(xs, ys, **options)
    decorate(xlabel="Standard normal") 

下面是分别对雄性和雌性企鹅的结果。

grouped = adelie.groupby("Sex")

weights_male = grouped.get_group("MALE")["Body Mass (g)"]
normal_probability_plot(weights_male, ls="--", label="Male")

weights_female = grouped.get_group("FEMALE")["Body Mass (g)"]
normal_probability_plot(weights_female, label="Female")

decorate(ylabel="Weight (g)") 

两组的正态概率图都接近一条直线，这表明重量的分布遵循正态分布。当我们把这两组放在一起时，它们的重量分布是两个具有不同均值的正态分布的混合——这样的混合并不总是由正态分布很好地建模。

现在，让我们考虑一些使正态分布对分析非常有用的数学特性。

正态分布

下面的类定义了一个表示正态分布的对象。它包含作为属性的参数mu和sigma2，它们分别代表分布的均值和方差。sigma2这个名字是一个提醒，方差是标准差的平方，通常表示为sigma。

class  Normal:
  """Represents a Normal distribution"""

    def  __init__(self, mu, sigma2):
  """Make a Normal object.

 mu: mean
 sigma2: variance
 """
        self.mu = mu
        self.sigma2 = sigma2

    def  __repr__(self):
  """Returns a string representation."""
        return f"Normal({self.mu}, {self.sigma2})"

    __str__ = __repr__ 

例如，我们将创建一个表示与雄性企鹅的重量具有相同均值和方差的正态分布的Normal对象。

m, s = weights_male.mean(), weights_male.std()
dist_male = Normal(m, s**2)
dist_male 

Normal(4043.4931506849316, 120278.25342465754) 

另一个与雌性企鹅的均值和方差相同的Normal对象。

m, s = weights_female.mean(), weights_female.std()
dist_female = Normal(m, s**2)
dist_female 

Normal(3368.8356164383563, 72565.63926940637) 

接下来，我们将向Normal类添加一个方法，该方法可以从正态分布中生成一个随机样本。为了向现有类添加方法，我们将使用一个 Jupyter 魔法命令add_method_to，该命令定义在thinkstats模块中。这个命令不是 Python 的一部分——它仅在 Jupyter 笔记本中工作。

%%add_method_to Normal

def sample(self, n):
    """Generate a random sample from this distribution."""
    sigma = np.sqrt(self.sigma2)
    return np.random.normal(self.mu, sigma, n) 

我们将使用sample来展示正态分布的第一个有用特性：如果你从两个正态分布中抽取值并将它们相加，总和的分布也是正态的。

例如，我们将从我们刚才创建的Normal对象中生成样本，将它们相加，并制作总和的正态概率图。

sample_sum = dist_male.sample(1000) + dist_female.sample(1000)
normal_probability_plot(sample_sum)

decorate(ylabel="Total weight (g)") 

正态概率图看起来像一条直线，这表明总和遵循正态分布。不仅如此——如果我们知道两个分布的参数，我们还可以计算总和分布的参数。以下方法展示了如何。

%%add_method_to Normal

def __add__(self, other):
    """Distribution of the sum of two normal distributions."""
    return Normal(self.mu + other.mu, self.sigma2 + other.sigma2) 

在总和的分布中，均值是均值的总和，方差是方差的和。现在我们已经定义了特殊方法__add__，我们可以使用+运算符来“相加”两个分布——也就是说，计算它们的总和的分布。

dist_sum = dist_male + dist_female
dist_sum 

Normal(7412.328767123288, 192843.8926940639) 

为了确认这个结果是正确的，我们将使用以下方法，它绘制了正态分布的解析 CDF。

%%add_method_to Normal

def plot_cdf(self, n_sigmas=3.5, **options):
    """Plot the CDF of this distribution."""
    mu, sigma = self.mu, np.sqrt(self.sigma2)
    low, high = mu - n_sigmas * sigma, mu + n_sigmas * sigma
    xs = np.linspace(low, high, 101)
    ys = norm.cdf(xs, mu, sigma)
    plt.plot(xs, ys, **options) 

这是结果，以及随机样本总和的经验 CDF。

dist_sum.plot_cdf(**model_options)
Cdf.from_seq(sample_sum).plot(label="sample")

decorate(xlabel="Total weight (g)", ylabel="CDF") 

看起来我们计算出的参数是正确的，这证实了我们可以通过相加均值和方差来“相加”两个正态分布。

作为推论，如果我们从一个正态分布中生成n个值并将它们相加，那么总和的分布也是一个正态分布。为了演示，我们将首先从男性体重的分布中生成 73 个值并将它们相加。下面的循环重复 1001 次，因此结果是总和分布的一个样本。

n = len(weights_male)
sample_sums_male = [dist_male.sample(n).sum() for i in range(1001)]
n 

73 

以下方法创建了一个Normal对象，它代表了总和的分布。为了计算参数，我们将均值和方差都乘以n。

%%add_method_to Normal

def sum(self, n):
    """Return the distribution of the sum of n values."""
    return Normal(n * self.mu, n * self.sigma2) 

这是n个重量的总和的分布。

dist_sums_male = dist_male.sum(n) 

这是它与随机样本的经验分布的比较。

dist_sums_male.plot_cdf(**model_options)
Cdf.from_seq(sample_sums_male).plot(label="sample")

decorate(xlabel="Total weights (g)", ylabel="CDF") 

解析分布与样本分布相匹配，这证实了sum方法是正确的。所以如果我们收集一个包含n个测量的样本，我们可以计算它们的总和的分布。

样本均值分布

如果我们可以计算样本总和的分布，我们也可以计算样本均值的分布。为了做到这一点，我们将使用正态分布的第三个性质：如果我们乘以或除以一个常数，结果仍然是一个正态分布。以下方法展示了我们如何计算乘积或商的分布的参数。

%%add_method_to Normal

def __mul__(self, factor):
    """Multiplies by a scalar."""
    return Normal(factor * self.mu, factor**2 * self.sigma2) 

%%add_method_to Normal

def __truediv__(self, factor):
    """Divides by a scalar."""
    return self * (1 / factor) 

为了计算乘积的分布，我们将均值乘以factor，将方差乘以factor的平方。我们可以使用这个性质来计算样本均值的分布。

dist_mean_male = dist_sums_male / n 

为了查看结果是否正确，我们还将计算随机样本的均值。

sample_means_male = np.array(sample_sums_male) / n 

将正常模型与样本均值的经验累积分布函数（CDF）进行比较。

dist_mean_male.plot_cdf(**model_options)
Cdf.from_seq(sample_means_male).plot(label="sample")

decorate(xlabel="Average weight (g)", ylabel="CDF") 

模型和模拟结果一致，这表明我们可以通过解析方法计算样本均值的分布——与重采样相比，这非常快。

现在我们知道了均值的抽样分布，我们可以用它来计算标准误差，这是抽样分布的标准差。

standard_error = np.sqrt(dist_mean_male.sigma2)
standard_error 

np.float64(40.591222045992765) 

这个结果建议了一个我们可以用来直接计算标准误差的捷径，而不需要计算抽样分布。在我们遵循的步骤中，我们首先将方差乘以n，然后除以n**2——净效果是将方差除以n，这意味着我们将标准差除以n的平方根。

因此，我们可以这样计算样本均值的标准误差。

standard_error = weights_male.std() / np.sqrt(n)
standard_error 

np.float64(40.59122204599277) 

现在，让我们考虑另一个我们可以用正态分布计算的结果，即差分的分布。

差分分布

将上一节中的步骤综合起来，以下是计算雌性企鹅体重样本均值分布的方法。

n = len(weights_female)
dist_mean_female = dist_female.sum(n) / n
dist_mean_female 

Normal(3368.835616438356, 994.0498530055667) 

现在我们有了雄性和雌性企鹅平均体重的抽样分布——让我们计算差分的分布。以下方法计算了两个正态分布值之间差分的分布。

%%add_method_to Normal

def __sub__(self, other):
    """Compute the distribution of a difference."""
    return Normal(self.mu - other.mu, self.sigma2 + other.sigma2) 

如你所料，差分的均值是均值的差。但正如你可能不会预料到的那样，差分的方差不是方差的差——它是总和！为了理解为什么，想象我们分两步进行减法：

• 如果我们对第二个分布取反，均值会被取反，但方差保持不变。

• 然后，如果我们加入第一个分布，总和的方差是方差的和。

如果这还不能说服你，让我们来测试它。以下是差分的解析分布。

dist_diff_means = dist_mean_male - dist_mean_female
dist_diff_means 

Normal(674.6575342465753, 2641.697160192656) 

这里是一个差分的随机样本。

sample_sums_female = [dist_female.sample(n).sum() for i in range(1001)]
sample_means_female = np.array(sample_sums_female) / n
sample_diff_means = sample_means_male - sample_means_female 

下图显示了随机样本的经验累积分布函数和正态分布的解析累积分布函数。

dist_diff_means.plot_cdf(**model_options)
Cdf.from_seq(sample_diff_means).plot(label="sample")

decorate(xlabel="Difference in average weight (g)", ylabel="CDF") 

它们一致，这证实了我们正确地找到了差分的分布。我们可以使用这个分布来计算体重差异的置信区间。我们将使用以下方法来计算逆累积分布函数。

%%add_method_to Normal

def ppf(self, xs):
    sigma = np.sqrt(self.sigma2)
    return norm.ppf(xs, self.mu, sigma) 

第 5 百分位数和第 95 百分位数形成一个 90%的置信区间。

ci90 = dist_diff_means.ppf([0.05, 0.95])
ci90 

array([590.1162635 , 759.19880499]) 

我们从随机样本中得到几乎相同的结果。

np.percentile(sample_diff_means, [5, 95]) 

array([578.86837628, 759.86657441]) 

解析方法比重采样更快，并且是确定的——也就是说，不是随机的。

然而，我们迄今为止所做的一切都是基于这样一个假设：测量的分布是正态的。这并不总是正确的——事实上，对于真实数据来说，它永远不是完全正确的。但即使测量的分布不是正态的，如果我们加起来很多测量值，它们的总和分布通常接近正态。这就是中心极限定理的力量。

中心极限定理

正如我们在前面的章节中看到的，如果我们加起来来自正态分布的值，总和的分布也是正态的。大多数其他分布没有这个特性——例如，如果我们加起来来自指数分布的值，总和的分布就不是指数分布。

但对于许多分布，如果我们生成n个值并将它们加起来，随着n的增加，总和的分布会收敛到正态分布。更具体地说，如果值的分布具有均值m和方差s2，总和的分布会收敛到具有均值n * m和方差n * s2的正态分布。

这个结论就是中心极限定理（CLT）。它是统计分析中最有用的工具之一，但同时也伴随着一些注意事项：

• 这些值必须来自同一个分布（尽管这个要求可以放宽）。

• 这些值必须独立抽取。如果它们是相关的，中心极限定理不适用（尽管如果相关性不是很强，它仍然可以工作）。

• 这些值必须来自具有有限均值和方差的分布。因此，中心极限定理不适用于一些长尾分布。

中心极限定理解释了正态分布在自然界中的普遍性。许多生物特征受到遗传和环境因素的累加效应的影响。我们测量的特征是大量小效应的总和，因此它们的分布倾向于正态分布。

为了了解中心极限定理是如何工作的，以及它在何时不起作用，让我们尝试一些实验，从指数分布开始。下面的循环从指数分布中生成样本，将它们加起来，并创建一个字典，将每个样本大小n映射到一个包含 1001 个总和的列表。

lam = 1
df_sample_expo = pd.DataFrame()
for n in [1, 10, 100]:
    df_sample_expo[n] = [np.sum(np.random.exponential(lam, n)) for _ in range(1001)] 

这里是每个总和列表的平均值。

df_sample_expo.mean() 

1        0.989885
10       9.825744
100    100.022555
dtype: float64 

这个分布的平均值是 1，所以如果我们加起来 10 个值，总和的平均值接近 10，如果我们加起来 100 个值，总和的平均值接近 100。

这个函数接受我们刚才制作的DataFrame，并为每个总和列表制作一个正常概率图。

def  normal_plot_samples(df_sample, ylabel=""):
  """Normal probability plots for samples of sums."""
    plt.figure(figsize=(6.8, 2.6))
    for i, n in enumerate(df_sample):
        plt.subplot(1, 3, i + 1)
        normal_probability_plot(df_sample[n])
        decorate(
            title="n=%d" % n,
            xticks=[],
            yticks=[],
            xlabel="Standard normal",
            ylabel=ylabel,
        ) 

下图显示了三个总和列表的正常概率图（normal_plot_samples的定义在本章的笔记本中）。

normal_plot_samples(df_sample_expo, ylabel="Sum of exponential values") 

当 n=1 时，总和的分布是指数分布，因此正态概率图不是一条直线。但是当 n=10 时，总和的分布近似为正态分布，当 n=100 时，几乎无法与正态分布区分。

对于比指数分布更倾斜的分布，总和的分布更快地收敛到正态分布——也就是说，对于较小的 n 值。对于更倾斜的分布，需要更长的时间。作为一个例子，让我们看看对数正态分布值的总和。

mu, sigma = 3.0, 1.0
df_sample_lognormal = pd.DataFrame()
for n in [1, 10, 100]:
    df_sample_lognormal[n] = [
        np.sum(np.random.lognormal(mu, sigma, n)) for _ in range(1001)
    ] 

下面是相同样本量范围的正态概率图。

normal_plot_samples(df_sample_lognormal, ylabel="Sum of lognormal values") 

当 n=1 时，正态模型不拟合分布，n=10 时也并不好多少。即使 n=100，分布的尾部也明显偏离模型。

对数正态分布的均值和方差是有限的，因此总和的分布最终会收敛到正态分布。但对于一些高度倾斜的分布，可能在任何实用的样本大小下都不会收敛。在某些情况下，甚至根本不会发生。

中心极限定理的局限性

帕累托分布比对数正态分布更加倾斜。根据参数的不同，一些帕累托分布没有有限的均值和方差——在这些情况下，中心极限定理不适用。

为了演示，我们将从参数 alpha=1 的帕累托分布中生成值，该分布具有无限的均值和方差。

alpha = 1.0
df_sample = pd.DataFrame()
for n in [1, 10, 100]:
    df_sample[n] = [np.sum(np.random.pareto(alpha, n)) for _ in range(1001)] 

下面是不同样本量下的正态概率图。

normal_plot_samples(df_sample, ylabel="Sum of Pareto values") 

即使 n=100，总和的分布也根本不像正态分布。

我还提到，如果值是相关的，中心极限定理不适用。为了测试这一点，我们将使用一个名为 generate_expo_correlated 的函数来生成指数分布的值，其中序列相关——即样本中连续元素之间的相关性——是给定的值，rho。此函数定义在本章的笔记本中。

def  generate_normal_correlated(n, rho):
  """Generates an array of correlated values from a standard normal dist."""
    xs = np.empty(n)
    xs[0] = np.random.normal(0, 1)

    sigma = np.sqrt(1 - rho**2)
    for i in range(1, n):
        xs[i] = rho * xs[i - 1] + np.random.normal(0, sigma)

    return xs 

给定一个来自正态分布的相关序列，以下函数生成一个来自指数分布的相关序列。

from  scipy.stats  import expon

def  generate_expo_correlated(n, rho):
  """Generates a sequence of correlated values from an exponential dist."""
    normal = generate_normal_correlated(n, rho)
    uniform = norm.cdf(normal)
    expo = expon.ppf(uniform)
    return expo 

它从一个相关的正态值序列开始，并使用正态 CDF 将它们转换为 0 到 1 之间的均匀分布的值序列。然后它使用指数逆 CDF 将它们转换为指数值序列。

以下循环创建一个 DataFrame，其中每列对应一个样本大小，每列有 1001 个总和。

rho = 0.8
df_sample = pd.DataFrame()
for n in [1, 10, 100]:
    df_sample[n] = [np.sum(generate_expo_correlated(n, rho)) for _ in range(1001)] 

下面是这些总和分布的正态概率图。

normal_plot_samples(df_sample, ylabel="Sum of correlated values") 

当 rho=0.8 时，相邻元素之间存在强烈的相关性，总和的分布缓慢收敛。如果序列中远距离元素之间也存在强烈的相关性，它可能根本不会收敛。

前一节展示了中心极限定理是如何工作的，本节展示了当它不起作用时会发生什么。现在让我们看看我们如何使用它。

应用中心极限定理

为了了解为什么中心极限定理是有用的，让我们回到第九章示例：测试第一胎和其他婴儿的平均怀孕长度之间的明显差异。我们再次使用 NSFG 数据——下载说明在本章的笔记本中。

以下单元格下载数据。

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") 

我们将使用 get_nsfg_groups 来读取数据并将其分为第一胎和其他婴儿。

from  nsfg  import get_nsfg_groups

live, firsts, others = get_nsfg_groups() 

正如我们所见，第一胎婴儿平均出生时间稍晚——明显的差异约为 0.078 周。

delta = firsts["prglngth"].mean() - others["prglngth"].mean()
delta 

np.float64(0.07803726677754952) 

为了看看这个差异是否可能是偶然发生的，我们将假设一个零假设，即怀孕长度的平均值和方差实际上对两组都是相同的，因此我们可以使用所有活产婴儿来估计它。

all_lengths = live["prglngth"]
m, s2 = all_lengths.mean(), all_lengths.var() 

怀孕长度的分布并不遵循正态分布——尽管如此，我们可以使用正态分布来近似均值的抽样分布。

以下函数接受一系列值并返回一个 Normal 对象，该对象表示从具有相同均值和方差的正态分布中抽取的给定大小 n 的样本均值的抽样分布。

def  sampling_dist_mean(data, n):
    mean, var = data.mean(), data.var()
    dist = Normal(mean, var)
    return dist.sum(n) / n 

这是根据零假设下第一胎出生的平均体重抽样分布的正常近似。

n1 = firsts["totalwgt_lb"].count()
dist_firsts = sampling_dist_mean(all_lengths, n1)
n1 

np.int64(4363) 

这是其他婴儿的抽样分布。

n2 = others["totalwgt_lb"].count()
dist_others = sampling_dist_mean(all_lengths, n2)
n2 

np.int64(4675) 

我们可以这样计算差异的抽样分布。

dist_diff = dist_firsts - dist_others
dist_diff 

Normal(0.0, 0.003235837567930557) 

均值为 0，这是有意义的，因为如果我们从同一分布中抽取两个样本，我们期望平均值的差异平均为 0。抽样分布的方差为 0.0032，这表明我们期望由于偶然性导致的差异变化有多大。

为了确认这个分布近似于抽样分布，我们还可以通过重采样来估计它。

sample_firsts = [np.random.choice(all_lengths, n1).mean() for i in range(1001)]
sample_others = [np.random.choice(all_lengths, n2).mean() for i in range(1001)]
sample_diffs = np.subtract(sample_firsts, sample_others) 

这是重采样差异的经验累积分布函数与正态模型相比。垂直虚线显示观察到的差异，正负都有。

dist_diff.plot_cdf(**model_options)
Cdf.from_seq(sample_diffs).plot(label="sample")
plt.axvline(delta, ls=":")
plt.axvline(-delta, ls=":")

decorate(xlabel="Difference in pregnancy length", ylabel="CDF") 

在这个例子中，样本量很大，测量的偏斜度适中，因此抽样分布很好地被正态分布所近似。因此，我们可以使用正态 CDF 来计算 p 值。以下方法计算正态分布的 CDF。

%%add_method_to Normal

def cdf(self, xs):
    sigma = np.sqrt(self.sigma2)
    return norm.cdf(xs, self.mu, sigma) 

这是根据零假设下差异达到 delta 的概率，即抽样分布右尾下的面积。

right = 1 - dist_diff.cdf(delta)
right 

np.float64(0.08505405315526993) 

以及这里负 -delta 差异的概率，这是左尾下的面积。

left = dist_diff.cdf(-delta)
left 

np.float64(0.08505405315526993) 

left 和 right 是相同的，因为正态分布是对称的。它们的和是 delta 差异发生的概率，无论是正还是负。

left + right 

np.float64(0.17010810631053985) 

得到的 p 值为 0.170，这与我们在第九章（chap09.html#section-diff-means）中通过重采样计算出的估计值一致。

我们计算这个 p 值的方法与独立样本 (t) 检验类似。SciPy 提供了一个名为 ttest_ind 的函数，它接受两个样本并计算它们均值差异的 p 值。

from  scipy.stats  import ttest_ind

result = ttest_ind(firsts["prglngth"], others["prglngth"])
result.pvalue 

np.float64(0.16755412639415004) 

当样本量较大时，(t) 检验的结果接近我们使用正态分布计算的结果。(t) 检验之所以称为 (t) 检验，是因为它基于 (t) 分布而不是正态分布。(t) 分布也用于测试相关性是否具有统计学意义，正如我们将在下一节中看到的。

相关性检验

在第九章（chap09.html#section-test-correlation）中，我们使用了排列检验来测试出生体重和母亲年龄之间的相关性，并发现它是统计显著的，p 值小于 0.001。

现在我们可以进行同样的分析。这种方法基于以下数学结果：如果我们从正态分布中生成两个大小为 n 的样本，计算皮尔逊相关系数 r，然后使用此函数转换相关性：

def  transform_correlation(r, n):
    return r * np.sqrt((n - 2) / (1 - r**2)) 

转换相关性遵循参数为 n-2 的 (t) 分布。为了了解这看起来像什么，我们将使用以下函数从标准正态分布生成不相关的样本。

def  generate_data(n):
  """Uncorrelated sequences from a standard normal."""
    xs = np.random.normal(0, 1, n)
    ys = np.random.normal(0, 1, n)
    return xs, ys 

以及以下函数来计算它们的相关性。

def  correlation(data):
    xs, ys = data
    return np.corrcoef(xs, ys)[0, 1] 

以下循环生成许多样本对，计算它们的关联性，并将结果放入列表中。

n = 100
rs = [correlation(generate_data(n)) for i in range(1001)] 

接下来，我们将计算转换相关性。

ts = transform_correlation(np.array(rs), n) 

为了检查这些 ts 是否遵循 (t) 分布，我们将使用以下函数，它创建一个表示 (t) 分布累积分布函数（CDF）的对象。

from  scipy.stats  import t as student_t

def  make_student_cdf(df):
  """Computes the CDF of a Student t distribution."""
    ts = np.linspace(-3, 3, 101)
    ps = student_t.cdf(ts, df=df)
    return Cdf(ps, ts) 

(t) 分布的参数称为 df，代表“自由度”。以下图显示了参数为 n-2 的 (t) 分布的 CDF 以及转换相关性的经验 CDF。

make_student_cdf(df=n - 2).plot(**model_options)

cdf_ts = Cdf.from_seq(ts)
cdf_ts.plot(label="random normals")

decorate(xlabel="Transformed correlation", ylabel="CDF") 

这表明如果我们从正态分布中抽取不相关的样本，它们的转换相关性将遵循 (t) 分布。

如果我们从其他分布中抽取样本，它们的转换相关性不会完全遵循 (t) 分布，但随着样本量的增加，它们会收敛到 (t) 分布。让我们看看这适用于母体年龄和出生体重的相关性。从活产婴儿的 DataFrame 中，我们将选择具有有效数据的行。

valid = live.dropna(subset=["agepreg", "totalwgt_lb"])
n = len(valid)
n 

9038 

实际的相关性约为 0.07。

data = valid["agepreg"].values, valid["totalwgt_lb"].values
r_actual = correlation(data)
r_actual 

np.float64(0.0688339703541091) 

如我们在 第九章 中所做的那样，我们可以通过重新排列样本来模拟零假设。

def  permute(data):
  """Shuffle the x values."""
    xs, ys = data
    new_xs = xs.copy()
    np.random.shuffle(new_xs)
    return new_xs, ys 

如果我们生成许多排列并计算它们的关联性，结果就是零假设下关联性分布的一个样本。

permuted_corrs = [correlation(permute(data)) for i in range(1001)] 

我们可以像这样计算转换后的相关性。

ts = transform_correlation(np.array(permuted_corrs), n) 

下图显示了 ts 的经验累积分布函数（CDF）以及参数为 n-2 的 t 分布的 CDF。

make_student_cdf(n - 2).plot(**model_options)
Cdf.from_seq(ts).plot(label="permuted data")

decorate(xlabel="Transformed correlation", ylabel="CDF") 

模型很好地拟合了经验分布，这意味着我们可以用它来计算观察到的相关性的 p 值。首先，我们将转换观察到的相关性。

t_actual = transform_correlation(r_actual, n) 

现在，我们可以使用 t 分布的 CDF 来计算在零假设下值为 t_actual 的概率。

right = 1 - student_t.cdf(t_actual, df=n - 2)
right 

np.float64(2.861466619208386e-11) 

我们还可以计算值为 -t_actual 的概率。

left = student_t.cdf(-t_actual, df=n - 2)
left 

np.float64(2.8614735536574016e-11) 

两者之和是相关系数为 r_actual 的概率，无论是正还是负。

left + right 

np.float64(5.722940172865787e-11) 

SciPy 提供了一个执行相同计算并返回观察到的相关性的 p 值的函数。

from  scipy.stats  import pearsonr

corr, p_value = pearsonr(*data)
p_value 

np.float64(5.7229471073151754e-11) 

结果几乎相同。

根据重采样结果，我们得出结论，p 值小于 0.001，但如果没有进行大量重采样，我们无法确定具体小多少。使用分析方法，我们可以快速计算小的 p 值。

然而，在实践中，这可能并不重要。一般来说，如果 p 值小于 0.001，我们可以得出结论，观察到的效应不太可能是由于偶然。通常，知道具体有多不可能并不重要。

卡方检验

在 第九章 中，我们基于这组观察到的结果测试了一个骰子是否是弯曲的。

from  empiricaldist  import Hist

qs = np.arange(1, 7)
freqs = [8, 9, 19, 5, 8, 11]
observed = Hist(freqs, qs)
observed.index.name = "outcome"
observed 

	频率
结果
---	---
1	8
2	9
3	19
4	5
5	8
6	11

首先，我们计算了每个结果的预期频率。

num_rolls = observed.sum()
outcomes = observed.qs
expected = Hist(num_rolls / 6, outcomes) 

然后，我们使用以下函数来计算卡方统计量。

def  chi_squared_stat(observed, expected):
    diffs = (observed - expected) ** 2
    ratios = diffs / expected
    return np.sum(ratios.values.flatten()) 

observed_chi2 = chi_squared_stat(observed, expected) 

卡方统计量在处理这类数据时被广泛使用，因为其在零假设下的抽样分布收敛到一个我们可以高效计算的分布——不是巧合，它被称为卡方分布。为了了解它的样子，我们将使用以下函数，该函数模拟掷一个公平的骰子。

def  simulate_dice(observed):
    n = np.sum(observed)
    rolls = np.random.choice(observed.qs, n, replace=True)
    hist = Hist.from_seq(rolls)
    return hist 

以下循环多次运行模拟并计算结果的卡方统计量。

simulated_chi_squared = [
    chi_squared_stat(simulate_dice(observed), expected) for i in range(1001)
]
cdf_simulated = Cdf.from_seq(simulated_chi_squared) 

为了检查结果是否遵循卡方分布，我们将使用以下函数，该函数计算参数为 df 的卡方分布的 CDF。

from  scipy.stats  import chi2 as chi2_dist

def  chi_squared_cdf(df):
  """Discrete approximation of the chi-squared CDF."""
    xs = np.linspace(0, 21, 101)
    ps = chi2_dist.cdf(xs, df=df)
    return Cdf(ps, xs) 

有 n 个可能的结果，模拟的卡方统计量应该遵循参数为 n-1 的卡方分布。

n = len(observed)
cdf_model = chi_squared_cdf(df=n - 1) 

这里是模拟的卡方统计量的经验 CDF 以及卡方分布的 CDF。

cdf_model.plot(**model_options)
cdf_simulated.plot(label="simulation")

decorate(xlabel="Chi-squared statistic", ylabel="CDF") 

模型很好地拟合了模拟结果，因此我们可以用它来计算在零假设下，值如observed_chi2那么大的概率。

p_value = 1 - chi2_dist.cdf(observed_chi2, df=n - 1)
p_value 

np.float64(0.04069938850404997) 

SciPy 提供了一个执行相同计算的功能。

from  scipy.stats  import chisquare

chi2_stat, p_value = chisquare(f_obs=observed, f_exp=expected) 

结果与我们所计算的 p 值相同。

p_value 

np.float64(0.040699388504049985) 

卡方统计量的优点是其零假设下的分布可以高效地计算。但在特定情况下，它可能不是最能量化观测值和预期值之间差异的统计量。

计算与分析

本书重点介绍计算方法，如重采样和排列。这些方法相对于分析具有几个优点：

• 它们更容易解释和理解。例如，在入门统计学课程中最难的主题之一是假设检验。许多学生并不真正理解 p 值是什么。我认为我们在第九章中采取的方法——模拟零假设并计算检验统计量——使基本思想更加清晰。

• 它们具有鲁棒性和多功能性。分析方法通常基于在实践中不成立的假设。计算方法需要的假设较少，并且更容易适应和扩展。

• 它们是可调试的。分析方法通常像是一个黑盒：你输入数字，它们就会输出结果。但是，很容易犯细微的错误，很难对结果有信心，如果结果不正确，也很难诊断问题。计算方法适合增量开发和测试，这有助于增强对结果的可信度。

但也存在缺点：

• 计算方法可能会很慢。

• 随机化方法，如重采样，每次产生的结果并不相同，这使得验证它们是否正确变得更加困难。

考虑到这些优缺点，我推荐以下流程：

1. 在探索过程中使用计算方法。如果你找到一个令人满意的答案，并且运行时间是可以接受的，那么你可以停止。

2. 如果运行时间不可接受，寻找优化机会。使用分析方法是几种优化方法之一。

3. 如果用分析方法替换计算方法合适，则将计算方法作为比较的基础，提供计算结果和分析结果之间的相互验证。

对于许多实际问题，计算方法的运行时间并不是问题，我们不需要进行到第一步之后。

术语表

• 正态概率图：一种比较观测值与正态分布分位数，以查看数据与正态分布的接近程度的图表。

• 独立样本 t 检验：一种计算两个独立组均值观测差异的 p 值的方法。

• t 分布：一种用于模拟在零假设下差异为 0 的均值差异的抽样分布以及转换后相关性的抽样分布。

• 卡方分布：一种用于模拟卡方统计量的抽样分布的分布。

• 卡方统计量：一个测试统计量，用于量化两个离散分布之间差异的大小。

练习

练习 14.1

在本章中，我们比较了雄性和雌性企鹅的体重，并计算了差异的置信区间。现在让我们对鳍长做同样的事情。观察到的差异大约是 4.6 毫米。

grouped = adelie.groupby("Sex")

lengths_male = grouped.get_group("MALE")["Flipper Length (mm)"]
lengths_female = grouped.get_group("FEMALE")["Flipper Length (mm)"]
observed_diff = lengths_male.mean() - lengths_female.mean()
observed_diff 

np.float64(4.616438356164366) 

使用sampling_dist_mean创建代表两组中平均鳍长抽样分布的Normal对象——注意两组的大小并不相同。然后计算差异的抽样分布和 90%的置信区间。

练习 14.2

使用 NSFG 数据，我们计算了婴儿出生体重与母亲年龄之间的相关性，并使用t分布来计算 p 值。现在让我们用出生体重和父亲的年龄来做同样的事情，这些年龄记录在hpagelb列中。

valid = live.dropna(subset=["hpagelb", "totalwgt_lb"])
n = len(valid)
n 

8933 

观察到的相关性约为 0.065。

data = valid["hpagelb"].values, valid["totalwgt_lb"].values
r_actual = correlation(data)
r_actual 

np.float64(0.06468629895432174) 

计算转换后的相关性，t_actual。使用t分布的累积分布函数（CDF）来计算 p 值——这个相关性是否具有统计学意义？使用 SciPy 函数pearsonr来检查你的结果。

练习 14.3

在第九章的一个练习中，我们考虑了特里弗斯-威拉德假说，该假说表明，对于许多哺乳动物，性别比例取决于“母体条件”——也就是说，像母亲的年龄、体型、健康状况和社会地位这样的因素。一些研究表明，在人类中也有这种效应，但结果不一。

作为例子，也是练习卡方检验的机会，让我们看看婴儿的性别和母亲的婚姻状况之间是否存在关系。本章节的笔记本中有指导说明，可以帮助你开始。

首先，我们将将男性和女性的母亲分开。

male = live.query("babysex == 1")
female = live.query("babysex == 2") 

现在我们将创建一个DataFrame，其中每一列代表一个组，每一行代表fmarital的每个值，如下所示：

1  married
2  widowed
3  divorces
4  separated
5  never married 

observed = pd.DataFrame()
observed["male"] = male["fmarital"].value_counts().sort_index()
observed["female"] = female["fmarital"].value_counts().sort_index()
observed 

	male	female
fmarital
---	---	---
1	2576	2559
2	56	54
3	568	572
4	355	330
5	1086	985

原假设是婚姻状况的分布在这两组中是相同的，因此我们可以使用整个数据集来计算它。

from  empiricaldist  import Pmf

pmf_fmarital = Pmf.from_seq(live["fmarital"])
pmf_fmarital 

	probs
fmarital
---	---
1	0.561653
2	0.012024
3	0.124617
4	0.075208
5	0.226498

为了计算期望值，我们将pmf_marital中的概率乘以每列中的总案例数。

expected = pd.DataFrame()
expected["male"] = pmf_fmarital * observed["male"].sum()
expected["female"] = pmf_fmarital * observed["female"].sum()
expected 

	male	female
fmarital
---	---	---
1	2606.630739	2527.437691
2	55.805641	54.110188
3	578.349366	560.778312
4	349.038916	338.434631
5	1051.175339	1019.239178

使用observed和expected来计算卡方统计量。然后使用卡方分布的 CDF 来计算 p 值。自由度应该是n-1，其中n是观察到的DataFrame中的值数。然后使用 SciPy 函数chisquared来计算卡方统计量和 p 值。提示：使用axis=None参数将整个DataFrame视为一个单独的测试，而不是每列一个测试。

这个测试是否为特里弗斯-威尔拉德假说提供了支持？

练习 14.4

我们在本章中用于分析组间差异的方法可以扩展到分析“差异的差异”，这是一种常见的实验设计。作为一个例子，我们将使用 2014 年一篇研究旨在减缓学生工程团队中性别刻板任务分配干预效果的论文中的数据。

Stein, L. A., Aragon, D., Moreno, D., & Goodman, J. (2014, October). 减缓学生工程团队中性别刻板任务分配干预的持久效应的证据。载于 2014 IEEE 前沿教育会议 (FIE) 论文集 (第 1-9 页)。IEEE。

可从 ieeexplore.ieee.org/document/7044435/ 获取。

在干预前后，学生回答了一份调查，要求他们在 7 点量表上对课堂项目每个方面的贡献进行评分。

在干预之前，男性学生在项目编程方面报告的分数高于女性学生：男性报告的平均分数为 3.57，标准误差为 0.28；女性报告的平均分数为 1.91，标准误差为 0.32。

干预之后，性别差距变小了：男性的平均分数为 3.44（SE 0.16）；女性的平均分数为 3.18（SE 0.16）。

1. 制作四个Normal对象来表示干预前后男性和女性学生的估计均值抽样分布。因为我们有估计均值的标准误差，所以我们不需要知道样本大小就可以得到抽样分布的参数。

2. 计算干预前后性别差距（即均值差异）的抽样分布。

3. 然后计算差异差异的抽样分布——即差距大小的变化。计算 95%置信区间和 p 值。

干预之后，性别差距的大小是否有所减小？

Think Stats: Python 中的探索性数据分析，第 3 版

版权所有 2024 艾伦·B·唐尼

代码许可：MIT 许可证

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