分析封锁对印度德里空气质量的影响#
你将做什么#
计算空气质量指数（AQI）并对其执行配对Student's t检验。
你将学到什么#
你将学习移动平均线的概念
你将学习如何计算空气质量指数（AQI）
你将学习如何执行配对Student's t检验并找到
t和p值你将学习如何解释这些值
你需要什么#
你的环境中安装了SciPy
对统计术语如总体、样本、均值、标准差等有基本理解
空气污染问题#
空气污染是我们面临的最突出的污染类型之一，它对我们的日常生活有直接影响。COVID-19大流行导致世界各地实行封锁；这提供了一个难得的机会来研究人类活动（或缺乏人类活动）对空气污染的影响。在本教程中，我们将研究德里（受空气污染影响最严重的城市之一）在2020年3月至6月封锁前后和封锁期间的空气质量。为此，我们首先根据收集到的污染物测量数据计算每小时的空气质量指数。接下来，我们将对这些指数进行抽样并执行配对Student's t检验。这将从统计学上向我们证明空气质量因封锁而改善，这与我们的直觉相符。
让我们首先将必要的库导入到我们的环境中。
import numpy as np
from numpy.random import default_rng
from scipy import stats
构建数据集#
我们将使用印度空气质量数据数据集的精简版。该数据集包含印度多个城市不同站点每小时和每天的空气质量数据和AQI（空气质量指数）。本教程提供的精简版包含德里从2019年5月31日至2020年6月30日的每小时污染物测量数据。它包含计算空气质量指数所需的标准污染物以及其他一些重要污染物的测量值：颗粒物（PM 2.5和PM 10）、二氧化氮（NO2）、氨（NH3）、二氧化硫（SO2）、一氧化碳（CO）、臭氧（O3）、氮氧化物（NOx）、一氧化氮（NO）、苯、甲苯和二甲苯。
让我们打印出前几行，以预览我们的数据集。
! head air-quality-data.csv
Datetime,PM2.5,PM10,NO2,NH3,SO2,CO,O3,NOx,NO,Benzene,Toluene,Xylene
2019-05-31 00:00:00,103.26,305.46,94.71,31.43,30.16,3.0,18.06,178.31,152.73,13.65,83.47,2.54
2019-05-31 01:00:00,104.47,309.14,74.66,34.08,27.02,1.69,18.65,106.5,79.98,11.35,76.79,2.91
2019-05-31 02:00:00,90.0,314.02,48.11,32.6,18.12,0.83,28.27,48.45,25.27,5.66,32.91,1.59
2019-05-31 03:00:00,78.01,356.14,45.45,30.21,16.78,0.79,27.47,44.22,21.5,3.6,21.41,0.78
2019-05-31 04:00:00,80.19,372.9,45.23,28.68,16.41,0.76,26.92,44.06,22.15,4.5,23.39,0.62
2019-05-31 05:00:00,83.59,389.97,39.49,27.71,17.42,0.76,28.71,39.33,21.04,3.25,23.59,0.56
2019-05-31 06:00:00,79.04,371.64,39.61,26.87,16.91,0.84,29.26,43.11,24.37,3.12,15.27,0.46
2019-05-31 07:00:00,77.32,361.88,42.63,27.26,17.86,0.96,27.07,48.22,28.81,3.32,14.42,0.41
2019-05-31 08:00:00,84.3,377.77,42.49,28.41,20.19,0.98,33.05,48.22,27.76,3.4,14.53,0.4
对于本教程的目的，我们只关注计算AQI所需的标准污染物，即PM 2.5、PM 10、NO2、NH3、SO2、CO和O3。因此，我们只会使用np.loadtxt导入这些特定列。然后，我们将切片并创建两个数据集：pollutants_A包含PM 2.5、PM 10、NO2、NH3和SO2，以及pollutants_B包含CO和O3。这两个数据集将略有不同地处理，我们稍后会看到。
pollutant_data = np.loadtxt("air-quality-data.csv", dtype=float, delimiter=",",
skiprows=1, usecols=range(1, 8))
pollutants_A = pollutant_data[:, 0:5]
pollutants_B = pollutant_data[:, 5:]
print(pollutants_A.shape)
print(pollutants_B.shape)
(9528, 5)
(9528, 2)
我们的数据集中可能包含缺失值，用 NaN 表示，所以让我们用 np.isfinite 快速检查一下。
np.all(np.isfinite(pollutant_data))
np.True_
至此，我们已成功导入数据并检查其完整性。让我们继续进行AQI计算！
计算空气质量指数#
我们将使用印度中央污染控制委员会采用的方法计算AQI。步骤总结如下：
收集标准污染物24小时平均浓度值；一氧化碳（CO）和臭氧（O3）为8小时平均浓度值。
使用以下公式计算这些污染物的子指数
\[ Ip = \dfrac{\text{IHi – ILo}}{\text{BPHi – BPLo}}\cdot{\text{Cp – BPLo}} + \text{ILo} \]其中，
Ip= 污染物p的子指数
Cp= 污染物p的平均浓度
BPHi= 浓度断点，即大于或等于Cp
BPLo= 浓度断点，即小于或等于Cp
IHi= 对应于BPHi的AQI值
ILo= 对应于BPLo的AQI值任意给定时间的最大子指数即为空气质量指数。
空气质量指数是借助断点范围计算的，如下图所示。
让我们创建两个数组来存储AQI范围和断点，以便我们稍后用于计算。
AQI = np.array([0, 51, 101, 201, 301, 401, 501])
breakpoints = {
'PM2.5': np.array([0, 31, 61, 91, 121, 251]),
'PM10': np.array([0, 51, 101, 251, 351, 431]),
'NO2': np.array([0, 41, 81, 181, 281, 401]),
'NH3': np.array([0, 201, 401, 801, 1201, 1801]),
'SO2': np.array([0, 41, 81, 381, 801, 1601]),
'CO': np.array([0, 1.1, 2.1, 10.1, 17.1, 35]),
'O3': np.array([0, 51, 101, 169, 209, 749])
}
移动平均线#
第一步，我们必须计算pollutants_A的24小时移动平均值，以及pollutants_B的8小时移动平均值。我们将使用np.cumsum和切片索引编写一个简单的函数moving_mean来实现这一目标。
为了确保这两个集合长度相同，我们将根据pollutants_A_24hr_avg的长度截断pollutants_B_8hr_avg。这也将确保我们在相同的时间段内拥有所有污染物的浓度。
def moving_mean(a, n):
ret = np.cumsum(a, dtype=float, axis=0)
ret[n:] = ret[n:] - ret[:-n]
return ret[n - 1:] / n
pollutants_A_24hr_avg = moving_mean(pollutants_A, 24)
pollutants_B_8hr_avg = moving_mean(pollutants_B, 8)[-(pollutants_A_24hr_avg.shape[0]):]
现在，我们可以使用np.concatenate连接这两个集合，形成一个包含所有平均浓度的单一数据集。请注意，我们必须按列连接数组，因此我们传递 axis=1 参数。
pollutants = np.concatenate((pollutants_A_24hr_avg, pollutants_B_8hr_avg), axis=1)
子指数#
每个污染物的子指数是根据AQI与标准断点范围之间的线性关系，按上述公式计算的
compute_indices 函数首先借助于我们上面创建的 AQI 和 breakpoints 数组，为输入浓度和污染物获取AQI类别的正确上限和下限以及断点浓度。然后，它将这些值输入到公式中以计算子指数。
def compute_indices(pol, con):
bp = breakpoints[pol]
if pol == 'CO':
inc = 0.1
else:
inc = 1
if bp[0] <= con < bp[1]:
Bl = bp[0]
Bh = bp[1] - inc
Ih = AQI[1] - inc
Il = AQI[0]
elif bp[1] <= con < bp[2]:
Bl = bp[1]
Bh = bp[2] - inc
Ih = AQI[2] - inc
Il = AQI[1]
elif bp[2] <= con < bp[3]:
Bl = bp[2]
Bh = bp[3] - inc
Ih = AQI[3] - inc
Il = AQI[2]
elif bp[3] <= con < bp[4]:
Bl = bp[3]
Bh = bp[4] - inc
Ih = AQI[4] - inc
Il = AQI[3]
elif bp[4] <= con < bp[5]:
Bl = bp[4]
Bh = bp[5] - inc
Ih = AQI[5] - inc
Il = AQI[4]
elif bp[5] <= con:
Bl = bp[5]
Bh = bp[5] + bp[4] - (2 * inc)
Ih = AQI[6]
Il = AQI[5]
else:
print("Concentration out of range!")
return ((Ih - Il) / (Bh - Bl)) * (con - Bl) + Il
我们将使用np.vectorize来利用向量化的概念。这简单地意味着我们不需要自己循环遍历污染物数组的每个元素。向量化是NumPy的关键优势之一。
vcompute_indices = np.vectorize(compute_indices)
通过为每个污染物调用我们的向量化函数 vcompute_indices，我们得到了子指数。为了获得原始形状的数组，我们使用 np.stack。
sub_indices = np.stack((vcompute_indices('PM2.5', pollutants[..., 0]),
vcompute_indices('PM10', pollutants[..., 1]),
vcompute_indices('NO2', pollutants[..., 2]),
vcompute_indices('NH3', pollutants[..., 3]),
vcompute_indices('SO2', pollutants[..., 4]),
vcompute_indices('CO', pollutants[..., 5]),
vcompute_indices('O3', pollutants[..., 6])), axis=1)
空气质量指数#
使用 np.max，我们找出每个时期的最大子指数，这就是我们的空气质量指数！
aqi_array = np.max(sub_indices, axis=1)
至此，我们得到了从2019年6月1日至2020年6月30日每小时的AQI。请注意，尽管我们最初使用的是5月31日的数据，但在移动平均步骤中我们对其进行了截断。
对AQI进行配对Student's t检验#
假设检验是一种描述性统计形式，用于帮助我们根据数据做出决策。从计算出的AQI数据中，我们想找出在封锁实施前后平均AQI是否存在统计学上的显著差异。我们将使用左尾配对Student's t检验来计算两个检验统计量——t 统计量 和 p 值。然后我们将这些值与相应的临界值进行比较以做出决策。
抽样#
我们现在将原始数据集中的 datetime 列导入到一个 datetime64 dtype 数组中。我们将使用此数组对AQI数组进行索引，并获取数据集的子集。
datetime = np.loadtxt("air-quality-data.csv", dtype='M8[h]', delimiter=",",
skiprows=1, usecols=(0, ))[-(pollutants_A_24hr_avg.shape[0]):]
由于德里从2020年3月24日开始全面封锁，因此封锁后子集的时间段为2020年3月24日至2020年6月30日。封锁前子集的时间长度与3月24日之前相同。
after_lock = aqi_array[np.where(datetime >= np.datetime64('2020-03-24T00'))]
before_lock = aqi_array[np.where(datetime <= np.datetime64('2020-03-21T00'))][-(after_lock.shape[0]):]
print(after_lock.shape)
print(before_lock.shape)
(2376,)
(2376,)
为了确保我们的样本 *近似地* 呈正态分布，我们选取样本量为 n = 30。before_sample 和 after_sample 是在全面封锁前后抽取的随机观测值集合。我们使用 random.Generator.choice 来生成样本。
rng = default_rng()
before_sample = rng.choice(before_lock, size=30, replace=False)
after_sample = rng.choice(after_lock, size=30, replace=False)
定义假设#
让我们假设封锁前后的样本均值之间没有显著差异。这将是零假设。备择假设将是均值之间*存在*显著差异，并且AQI*有所改善*。数学上表示为：
\(H_{0}: \mu_\text{after-before} = 0\)
\(H_{a}: \mu_\text{after-before} < 0\)
计算检验统计量#
我们将使用 t 统计量来评估我们的假设，甚至从中计算 p 值。t 统计量的公式是
其中，
\(\mu_\text{after-before}\) = 样本均值差
\(\sigma^{2}\) = 均值差的方差
\(n\) = 样本量
def t_test(x, y):
diff = y - x
var = np.var(diff, ddof=1)
num = np.mean(diff)
denom = np.sqrt(var / len(x))
return np.divide(num, denom)
t_value = t_test(before_sample, after_sample)
对于 p 值，我们将使用 SciPy 的 stats.distributions.t.cdf() 函数。它接受两个参数：t 统计量 和自由度 (dof)。dof 的公式是 n - 1。
dof = len(before_sample) - 1
p_value = stats.distributions.t.cdf(t_value, dof)
print("The t value is {} and the p value is {}.".format(t_value, p_value))
The t value is -7.4504287940960126 and the p value is 1.640106293121364e-08.
t 和 p 值意味着什么？#
我们现在将计算出的检验统计量与临界检验统计量进行比较。临界 t 值通过查阅 t分布表 获得。
从上表中可以看出，当自由度为29且置信水平为95%时，临界值为1.699。由于我们使用的是左尾检验，因此我们的临界值为-1.699。显然，计算出的 t 值小于临界值，因此我们可以安全地拒绝零假设。
临界 p 值，由 \(\alpha\) 表示，通常选择为0.05，对应于95%的置信水平。如果计算出的 p 值 小于 \(\alpha\)，则可以安全地拒绝零假设。显然，我们的 p 值 远小于 \(\alpha\)，因此我们可以拒绝零假设。
请注意，这并不意味着我们可以接受备择假设。它只告诉我们没有足够的证据拒绝 \(H_{a}\)。换句话说，我们未能拒绝备择假设，因此它 *可能* 是真的。
在实践中…#
对于时间序列数据分析，更推荐使用 pandas 库。
SciPy 的 stats 模块提供了 stats.ttest_rel 函数，可用于获取
t 统计量和p 值。在实际生活中，数据通常不呈正态分布。对于此类非正态数据，存在像 Wilcoxon 检验 这样的测试。
进一步阅读#
你可以根据给定数据的特点选择多种统计检验。在 统计数据分布入门 中阅读更多相关内容。
你可以根据自己的需求采用 Student's t检验 的各种版本。