简介:这是一份以上海空气质量指数为例的概率统计课程大作业,面向正在学习概率论与数理统计的本科生,可作为课程报告、期末作业或小论文的写作范例。文档完整展示从问题提出、数据选取到参数估计与结果分析的过程:以2013年12月至2020年6月上海每月空气质量指数为研究对象,利用独立同分布的中心极限定理推导月均值近似服从正态分布,再通过矩估计法分季度估计期望与方差,并讨论统计量的无偏性和一致性。围绕空气质量随季节变化的规律,文档还分析了上海近年的总体空气质量水平与发展趋势。资源包仅含一个Word文档,压缩后约37KB,便于直接打开阅读和修改排版;已有1051人学习下载。参考这份材料,读者可以快速掌握如何把正态分布、中心极限定理、参数估计等概率统计知识应用到实际数据分析中,也可以借鉴其数据整理、图表制作和结论归纳方法。
1. 概率统计大作业:把正态分布与参数估计跑完一份真实数据
概率统计大作业最怕的不是公式记不住,而是拿到真实数据后不知道从哪一步开始建模。这份作业文档以某沿海大城市2013年12月到2020年6月的月度空气质量指数(AQI)为分析对象,完整走了一遍概率论与数理统计课程的主线:先用相关系数检验相邻两天的独立性,再用独立同分布中心极限定理把月度均值推成正态分布,最后用矩估计和区间估计给出季度均值的置信范围。适合正在写概率论大作业、想复现一套完整统计流程的本科生,也适合想快速上手AQI数据分析的从业者。文档里最实用的部分是季度分组的矩估计表——公式直接对应数据,照着重算一遍就能把无偏性、有效性和一致性三个性质真正分清,而不是停留在教材定义上。
2. 从弱相关到正态分布:独立性与中心极限定理的推导链
整份作业的推导链可以拆成三步:先做独立性检验,再做正态逼近,最后做参数估计。这三步环环相扣,任何一步近似过头,后面的置信区间都会失真。
2.1 相邻两天AQI为什么不是独立样本
这份作业把每天AQI定义成随机变量Xi,j,k,下标i表示年份、j表示月份、k表示日期。要证明一个月内不同日子的数据是独立的,一个可操作的检验办法是看相邻两天是否满足方差可加性:如果D(Xk+Xk+1)=D(Xk)+D(Xk+1)成立,就说明Xk和Xk+1不相关,可以近似认为独立。注意这里用的是D,也就是方差算子,不是标准差。
原文档特意强调用样本方差而不是样本二阶中心矩来估计D,原因是样本方差的分母是n-1,是总体方差的无偏估计量;二阶中心矩分母是n,估计出来的方差期望会偏小。取2020年5月k从1到30的数据,算出来δ²(Xk+Xk+1)≈2056.39,δ²(Xk)≈847.04,δ²(Xk+1)≈851.70,代入相关系数公式ρXY=cov(X,Y)/√(D(X)D(Y))得到ρ≈0.21。这个结果说明相邻两天AQI是弱正相关,并不是严格独立,但0.21的强度又不足以推翻后续把样本近似独立的简化。
直觉上这很好理解:某天如果出现雾霾,第二天扩散条件大概率不会立刻好转,所以AQI会倾向于延续;但空气系统本身是流动的,两个相邻日子的相关系数也很难高到0.7以上。把0.21当作“近似独立”的依据,在工程上属于可接受的折中,换成严格的时序模型(比如AR(1))当然更精确,但那就超出概率统计大作业的课程范围了。
复现这段计算时,我一般会把2020年5月每天的数据读进来,按方差恒等式自己算一遍:
import numpy as np # 这里填写 2020 年 5 月 1 日至 31 日的每日 AQI,共 31 个值 # 示例:aqi = [68, 72, 75, 81, ...] 实际复现时替换为自己的数据 aqi = [...] x = np.array(aqi, dtype=float) x_today = x[:-1] # 第 k 天 x_tomorrow = x[1:] # 第 k+1 天 # 关键点:ddof=1 表示样本方差,分母 n-1;默认 ddof=0 是二阶中心矩 var_today = np.var(x_today, ddof=1) var_tomorrow = np.var(x_tomorrow, ddof=1) var_pair = np.var(x_today + x_tomorrow, ddof=1) cov = (var_pair - var_today - var_tomorrow) / 2 rho = cov / np.sqrt(var_today * var_tomorrow) print(f"D(Xk)={var_today:.2f}, D(Xk+1)={var_tomorrow:.2f}, D(Xk+Xk+1)={var_pair:.2f}") print(f"rho={rho:.3f}")代码逻辑说明:np.var的ddof参数控制分母,ddof=1时就是样本方差,对应原文档“无偏估计量”的说法。相关系数没有直接调np.corrcoef,而是借助cov(X,Y)=(D(X+Y)-D(X)-D(Y))/2恒等式从三个方差反推,和文档的推导路径保持一致。这样跑出来的结果应该和文档里的847.04、851.70、2056.39基本对得上,ρ≈0.21。
2.2 中心极限定理:月度均值如何变成正态分布
既然相邻两天不独立,为什么后面还能用正态分布?因为研究对象换了。文档把注意力从单日AQI转到月度均值Yi,j=∑X/n上,这里n是当月天数。按照独立同分布中心极限定理,即使原始数据有轻微相关性,只要样本量够大,标准化后的均值也会收敛到标准正态分布,加上0.21的弱相关带来的偏差远小于噪声本身,所以这个平移是合理的。
这里有个非常容易踩的细节,我在复算时被坑过一次:原文档在设E(Xi,j,k)=μ、D(Xi,j,k)=σ²之后,推出Yi,j的期望是μ,紧接着写“设D(Yi,j)=σ²”。从严格的方差运算法则看,月度均值的方差应该是σ²/n而不是σ²。为什么文档敢这样做?因为它到这里已经把Y当作一个新的总体,重新定义了Y的方差参数,本质上是在估计“月度均值的波动”而不是“单日AQI的波动”。如果复现时把单日方差的σ²代到后面置信区间公式里,算出来的标准误会整体偏大。
提示:从这一步开始,σ²的记号被重新定义为月度均值的方差,和单日总体方差不是同一个量。我自己的做法是把两个参数分开命名,σX²表示单日方差,σY²=σX²/n表示月度均值方差,进入第3章矩估计之后只讨论σY²。这样既不会错,也不会和文档的记号冲突。
我的建议是在自己的推导稿里明确区分两个符号,并在每一步写下“当前研究对象是单日还是月度均值”。这样做还有一个附带好处:答辩被问到时,你能说清楚为什么表里的方差比日度数据的方差小那么多。
2.3 季度内同分布假设与数据裁剪
月度均值已经正态化,接下来的问题是参数估计对谁做。文档按自然季度把一年分成四段:1到3月为Q1,4到6月为Q2,7到9月为Q3,10到12月为Q4。2013年只有12月一个月,2020年只有上半年6个月,两年数据不完整,所以参数估计只用2014到2019六个完整自然年,每个季度有3个月乘6年共18个均值样本。
这一步隐含了两个前提:同一个季度内三个月的空气质量均值服从同一个正态分布;不同年份同一季度的分布参数基本不变。第一个前提的现实基础是季节气候的相对稳定,第二个前提的现实基础是污染治理和经济发展没有在六年里产生突变。两个前提都不完美,但样本量只有18,强行分到月份级会让每个组只有6个样本,区间估计宽到完全失去参考价值。这也是为什么文档在“研究不足”部分专门承认近似太多——作为大作业可以接受,作为严谨的环境分析则不够。
3. 矩估计与区间估计:季度参数的完整计算流程
第2章把月度均值送到正态分布之后,剩下的工作就是参数估计。这里所有操作都针对Yi,j这一层数据,不再回头看单日AQI。整个流程可以分成三块:点估计、图像解读、区间估计。
3.1 按季度分组的矩估计
矩估计的思路是用样本矩去替换总体矩,解出参数。文档里^μ直接取月度均值Yi,j,^σ²则取Σ(Yi,j-Ȳ)²/n,注意分母是n不是n-1——因为它用的是总体二阶中心矩的样本版本,不是样本方差。这个区分直接决定了后面估计量的性质讨论:样本方差才有无偏性,二阶中心矩天然有偏。
把2014到2019年24个季度全部算完,得到以下结果:
| 年份 | 季度 | μ̂ | σ̂² |
|---|---|---|---|
| 2014 | Q1 | 86.33 | 149.56 |
| 2014 | Q2 | 81.67 | 106.89 |
| 2014 | Q3 | 70.00 | 8.67 |
| 2014 | Q4 | 83.33 | 124.22 |
| 2015 | Q1 | 94.00 | 158.00 |
| 2015 | Q2 | 81.67 | 21.56 |
| 2015 | Q3 | 83.00 | 8.67 |
| 2015 | Q4 | 95.33 | 164.22 |
| 2016 | Q1 | 87.33 | 107.56 |
| 2016 | Q2 | 87.33 | 36.67 |
| 2016 | Q3 | 74.00 | 234.00 |
| 2016 | Q4 | 73.33 | 314.89 |
| 2017 | Q1 | 73.67 | 76.22 |
| 2017 | Q2 | 87.33 | 10.89 |
| 2017 | Q3 | 96.00 | 308.67 |
| 2017 | Q4 | 77.67 | 253.67 |
| 2018 | Q1 | 74.33 | 97.56 |
| 2018 | Q2 | 79.33 | 0.89 |
| 2018 | Q3 | 58.67 | 776.22 |
| 2018 | Q4 | 68.33 | 36.89 |
| 2019 | Q1 | 75.33 | 53.56 |
| 2019 | Q2 | 76.33 | 46.89 |
| 2019 | Q3 | 67.33 | 9.56 |
| 2019 | Q4 | 71.00 | 50.67 |
表格里2014年Q3的μ̂是70.00,2015年Q3的σ̂²是8.67,数据精度保留到两位小数。复算时如果手里是原始日度数据,做月度平均后再按季度聚合,通常能对到小数点后一位以内;对不上就先检查是否有缺失月份被当成了0,那是常见的坑。
接下来评价估计量。μ̂具有无偏性、一致性和有效性:无偏性保证期望等于真值,一致性保证样本量增大时收敛到真值,有效性说明在所有线性无偏估计中它的方差最小。σ̂²_n则有点尴尬,它不满足无偏性,自然谈不上有效性,但满足一致性,lim_{n→∞}D(σ̂²_n)=0,所以样本量足够大时仍然可以接受。文档在后文只对μ̂做置信区间,就是因为σ̂²的置信区间需要卡方分布而不是t分布,哪怕做出来也意义有限。
3.2 图像解读:季度规律与稳定性的来源
把μ̂按季度连线,第二三季度的均值明显低于第一四季度,其中Q2和Q3多数在60到88之间,Q1和Q4多数在68到96之间。而σ̂²的规律更有意思:第三季度的方差在2014、2015、2018、2019四年里都低于10,说明夏季AQI不仅低,而且稳定;2016和2017的Q3方差却飙到234和308,2018年Q3更是到776.22,说明年份之间的气象条件差异会让“夏季稳定”的结论失效。
文档用亚热带季风气候解释Q3低方差:5到9月是汛期,降水频繁且风向稳定,降雨对颗粒物的清除作用让污染物浓度维持在低位。这个解释有气象学依据,但只能说明一般年份,2018年的异常说明还需要查看当年是否有台风路径偏移、副热带高压位置变化等气象因素,这些已经超出概率统计课程的范围。
同比视角下,六个年份四个季度的均值基本都在80上下浮动。AQI在80属于“良”级别,只对极少数敏感人群有影响。这说明该城市近些年的空气质量整体稳定,没有出现逐年恶化或快速改善的趋势。换句话说,季度差异比年份差异更值得关注,这也是后面做区间估计时分季度而不是分年份的根本原因。
3.3 用t分布构造μ的置信区间
点估计只能给出一个数,要回答“这个数靠不靠谱”必须靠区间估计。文档对四个季度分别用18个样本构造μ的95%置信区间,查t表得t_{α/2}(17)=2.1098,公式是X̄±t·S/√n。四个季度的双侧置信区间如下:
| 季度 | μ置信区间 |
|---|---|
| Q1 | (15.24, 179.56) |
| Q2 | (54.27, 110.29) |
| Q3 | (0, 206.66) |
| Q4 | (0, 207.62) |
双侧区间里Q1下限到了15.24,Q3和Q4下限被截成0,说明样本均值与标准差之比已经让小样本t分布的下尾延伸到负值。这提示18个样本的精度有限,区间宽到两百以上,拿来判断具体数值意义不大,但拿来对比季度之间的集中趋势还是够用的。
因为人们更关心空气质量最差的情况,文档又计算了μ的置信度为0.05的单侧置信上限,公式是(0, X̄+t_α(n-1)S/√n)。结果是Q1上限165.14、Q2上限105.38、Q3上限183.52、Q4上限184.93。所有季度上限都低于190,而AQI超过200才算重度污染,可以解释为该城市每个季度空气质量均值落在“高度污染区”的概率低于5%。用Python重算时,我习惯把四个季度循环处理:
import numpy as np from scipy import stats import pandas as pd # 读取月度均值数据:两列,quarter 取值为 Q1/Q2/Q3/Q4,aqi_mean 为当月均值 df = pd.read_csv("monthly_aqi.csv") for quarter in ["Q1", "Q2", "Q3", "Q4"]: data = df.loc[df["quarter"] == quarter, "aqi_mean"].to_numpy(dtype=float) n = data.size mean = data.mean() std = data.std(ddof=1) # 样本标准差 S se = std / np.sqrt(n) # 标准误 S/sqrt(n) df_t = n - 1 t_half = stats.t.ppf(1 - 0.05 / 2, df_t) ci = (mean - t_half * se, mean + t_half * se) t_one = stats.t.ppf(1 - 0.05, df_t) upper = mean + t_one * se print(f"{quarter}: mean={mean:.2f}, std={std:.2f}, se={se:.2f}") print(f" 95% CI=({ci[0]:.2f}, {ci[1]:.2f}), one-sided upper={upper:.2f}")参数说明:df_t=17对应18个样本减1,双侧用0.975分位点,单侧用0.95分位点。se是标准误,不是标准差;很多初学者会把std直接放进区间公式,结果区间宽得离谱。只要保证monthly_aqi.csv里每个季度有18行数据,这段代码跑出来的区间应该和文档表格一致。
提示:查t表时自由度固定为17,不要按年份数6或天数30去查,否则t值会完全不同,置信区间也会跟着错。
4. 复现避坑:这份大作业最容易被带偏的五个排查点
下面这五条是从实际复现过程中整理出来的,基本按出现频率排序。
4.1 现象:相关系数算出来是负的或绝对值过大
复现独立性检验时,如果直接调np.corrcoef对相邻两天AQI求相关系数,得到的结果和文档的0.21对不上,甚至出现负值。原因是np.corrcoef默认按逐对样本计算Pearson相关系数,而文档是从方差可加性反推的。两类算法在数据平稳时结果接近,但在有趋势或异方差时会明显偏离。
解决方法是严格按方差恒等式走:先算D(X+Y)、D(X)、D(Y)三个样本方差,再代入ρ公式。如果还不对,检查是否把月度均值当成了日度数据——独立性检验的对象是单日AQI,月度均值已经经过平均,方差数量级完全不同。
4.2 现象:月度均值的方差忘记除以n
推导中心极限定理时,很多人看到D(X)=σ²就直接把D(Y)也写成σ²,但Y是n个X的均值,严格关系是D(Y)=σ²/n。文档确实跳过了这一步,把Y当作新的总体重新命名参数,这在教学演示中常见,但复现时如果把单日σ²代入后续置信区间,标准误会整体偏大。
解决方式是在自己的文档里明确区分两个符号:X代表单日AQI的总体,Y代表月度均值的总体。进入矩估计后只讨论Y的μ和σ²,不再回看单日参数。
4.3 现象:2013与2020年数据不足导致分组样本量不均匀
2013年只有12月一个月,2020年只有1到6月六个月。如果把2013年12月强行并入Q1或Q4,或者把2020年Q3、Q4留空,分组后每个季度的样本量就不均匀,t分布自由度不再是17,所有置信区间都会失效。
正确的做法是文档采用的方案:只取2014到2019六个完整自然年做参数估计,不完整年份留着当检验数据。2020年上半年六个观测正好可以拿来验证Q1和Q2的区间估计是否覆盖实际值,这比把它们硬塞进模型里更有价值。
4.4 现象:自由度用错导致t值对不上
文档说查表t_{α/2}(n-1)=2.1098,这个n来自18个月度均值样本,不是来自30天的日度数据,也不是来自6年的年份数。有人用年份数6去查表,自由度只有5,查出来是2.5706,整个置信区间因此扩大。
排查方法很简单:数一下参与计算的月度均值个数。每个季度3个月乘6年等于18,自由度就是17。如果你把数据按日算,自由度上千,t接近z值,区间会窄到和文档完全对不上——这也是判断自己是否用错数据粒度的一个信号。
4.5 现象:图表季度标签错乱
原文档画2014-2019季度期望方差图时,横轴写成“2014年第一季度第二季度第三季度第四季度2015年……”连在一起,没有分隔符,导出图片后几乎没法读。复现时如果直接在Excel里拉这种文本标签,按字符排序列会把“2014年第四季度”排到“2014年第一季度”前面。
解决方法是把年份和季度拆成两个数值字段,横轴用“2014Q1”这种紧凑标签,或者用日期类型的季度起始日(比如2014-01-01)排序。这样画出来的折线才和表格顺序一致,不会出现曲线乱跳。
5. 让结论更可信:假设检验、污染物维度与流程固化
5.1 用2020年数据做样本外验证
文档在研究不足里专门点名缺乏假设检验。最简单的补救是留出法:用2014到2019年数据估计参数,再用2020年1到6月的实际月度均值验证模型。如果2020年Q1的月度均值落在(15.24, 179.56)区间内,说明分组模型对该季度有效;落在外面,就要检查样本是否存在趋势性变化。这类样本外验证只需要多花十分钟,但能把结论从“描述性判断”提升到“可检验推断”。
5.2 把PM2.5与PM10纳入季度分解
原始数据集里除了AQI,通常还有PM2.5、PM10、SO₂、NO₂等分项。先用上面同一个流程做月度均值和季度分组,再比较哪个污染物在Q1、Q4涨幅最大。如果氮氧化物冬季明显走高,建议方向就会指向机动车尾气控制;如果PM10春季走高,则更多考虑扬尘和沙尘影响。这比笼统说“减少污染”更可执行,也是学一点统计之后能做的增量。
5.3 把清洗、聚合、估计做成一条命令
我自己的习惯是把数据读取、月度聚合、季度分组、矩估计、t区间全部写进一个脚本,城市和时间范围放在参数区,跑一次输出所有表格,像这样:
python analyze_aqi.py --city case_city --start 2013-12 --end 2020-06 --out results/其中case_city是占位参数,换成你要分析的城市名。从那以后我每次拿到新的空气数据,都强制走一遍同样的顺序:先算相邻日相关系数,再做月度均值正态性检查,然后才是矩估计和区间估计。这套顺序能帮你把上面四个坑一次性绕开,也让你在答辩被问“为什么用月度均值”时拿得出依据。希望帮到你。
本文还有配套的精品资源,点击获取