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

网格显示印度门上方烟雾缭绕，下方空气清澈

你会做什么#

计算空气质量指数 (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导入这些特定列。然后，我们将切片并创建两个集合：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)

我们的数据集可能包含缺失值，用表示，因此让我们使用np.isfiniteNaN进行快速检查。

np.all(np.isfinite(pollutant_data))

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= 对应的AQI值BPHi
ILo= 对应的AQI值BPLo
任何给定时间的最大分项指数是空气质量指数。

空气质量指数是借助断点范围计算的，如下图所示。

断点范围图表

让我们创建两个数组来存储 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])
}

移动平均线#

第一步，我们必须计算24 小时窗口和8 小时窗口的移动平均值。我们将使用np.cumsum和切片索引编写一个简单的函数来实现此目的。pollutants_Apollutants_Bmoving_mean

为了确保两个集合的长度相同，我们将pollutants_B_8hr_avg根据的长度截断pollutants_A_24hr_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与标准断点范围的线性关系计算，公式如下：

\[ Ip = \dfrac{\text{IHi – ILo}}{\text{BPHi – BPLo}}\cdot{\text{Cp – BPLo}} + \text{ILo} \]

该函数首先借助我们上面创建compute_indices的数组获取输入浓度和污染物的 AQI 类别的正确上限和下限以及断点浓度。然后，它将这些值输入公式来计算分项指数。AQIbreakpoints

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 的配对学生 t 检验#

假设检验是描述性统计的一种形式，用于帮助我们利用数据做出决策。从计算的 AQI 数据中，我们想了解封锁前后平均 AQI 是否存在统计上的显着差异。我们将使用左尾配对学生 t 检验来计算两个检验统计量- 和。然后我们将这些与相应的临界值进行比较以做出决定。t statisticp value

正态分布图显示单尾测试中的拒绝面积（左尾）

采样#

现在，我们将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,)

为了确保我们的样本近似正态分布，我们抽取大小为的样本。和是在完全封锁之前和之后绘制的一组随机观察值。我们使用random.Generator.choice来生成样本。n = 30before_sampleafter_sample

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 valuet

\[ t = \frac{\mu_\text{after-before}}{\sqrt{\sigma^{2}/n}} \]

在哪里，

\(\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 statisticdofdofn - 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.244515326256721 and the p value is 2.8076806429430512e-08.

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

我们现在将计算的测试统计数据与关键测试统计数据进行比较。临界值t是通过查t分布表计算出来的。

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

dof从上表中可以看出，置信水平为 95% 时，29 的临界值为 1.699 。由于我们使用左尾检验，因此我们的临界值为-1.699。显然，计算t值小于临界值，因此我们可以安全地拒绝原假设。

临界p值，表示为\(\alpha\)，通常选择 0.05，对应于 95% 的置信水平。如果计算p值小于\(\alpha\)，那么可以安全地拒绝原假设。显然，我们的p价值远低于\(\alpha\)，所以我们可以拒绝原假设。

请注意，这并不意味着我们可以接受备择假设。它只是告诉我们没有足够的证据来拒绝\(H_{a}\)。换句话说，我们无法拒绝替代假设，因此它可能是正确的。

在实践中… ＃

pandas库更适合用于时间序列数据分析。
SciPy stats 模块提供stats.ttest_rel函数，可用于获取和。t statisticp value
在现实生活中，数据通常不呈正态分布。有针对此类非正态数据的检验，例如Wilcoxon 检验。

进一步阅读#

您可以根据给定数据的特征选择许多统计检验。有关它们的更多信息，请参阅统计数据分布的简要介绍。
学生 t 检验有多种版本，您可以根据需要采用。

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

内容