分析封锁对印度德里空气质量的影响#

显示雾霾中的印度门在上，下方是清晰空气的网格

你要做什么#

计算空气质量指数 (AQI) 并对它们进行配对的学生 t 检验。

你将学到什么#

你将学习移动平均线的概念
你将学习如何计算空气质量指数（AQI）
你将学习如何执行配对的学生 t 检验并找到 t 和 p 值
你将学习如何解释这些值

你需要什么#

SciPy 已安装在你的环境中
对统计学术语的基本理解，如总体、样本、均值、标准差等。

空气污染的问题#

空气污染是我们面临的最突出的污染类型之一，它对我们的日常生活有直接影响。COVID-19 大流行导致世界不同地区的封锁，提供了一个难得的机会来研究人类活动（或缺乏活动）对空气污染的影响。在本教程中，我们将研究德里在 2020 年 3 月至 6 月封锁前后的空气质量，德里是受空气污染影响最严重的城市之一。为此，我们将首先根据收集到的污染物测量值计算每小时的空气质量指数。接下来，我们将对这些指数进行采样，并对它们进行配对样本 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。然后，我们将切片并创建两个集合：包含PM 2.5、PM 10、NO2、NH3和SO2的pollutants_A，以及包含CO和O3的pollutants_B。这两个集合将稍有不同的处理方式，我们将在后面看到。

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 = 污染物 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小时的移动平均。我们将编写一个简单的函数 moving_mean 使用 np.cumsum 和切片索引来实现这一点。

为了确保两组数据的长度相同，我们将根据 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日的数据开始的，但在移动平均步骤中我们截断了那些数据。

配对学生 t 检验在 AQI 上#

假设检验是一种描述性统计的形式，用于帮助我们根据数据做出决策。从计算的AQI数据中，我们想找出在实施封锁之前和之后，平均AQI是否有统计上的显著差异。我们将使用左尾的配对学生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 -9.523339365695536 and the p value is 9.875512270224618e-11.

`t` 和 `p` 值是什么意思？#

我们现在将比较计算的检验统计量与临界检验统计量。临界 t 值是通过查找 t 分布表计算的。

在不同置信水平下选择的t值表。在95%置信水平下29自由度的t值用黄色方块突出显示

从上表中，对于29个 dof 在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 检验。

进一步阅读#

根据给定数据的特点，你可以选择多种统计测试。更多信息请阅读 A Gentle Introduction to Statistical Data Distributions。
根据您的需求，您可以选择不同版本的 Student’s t-test。

分析封锁对印度德里空气质量的影响

目录