《数据分析基础——基于python的实现》笔记

统计基础

中心极限定理（Central Limit Theorem）

不知道为啥我看到的中心极限定理有两个版本的表述
（后来发现确实是有两个版本）

第一个版本说：某城市的工资分布是个很奇怪的分布
但如果对该城市进行抽样，每次抽20个人求平均值，抽100次，那么这100个平均值的分布就会是正态分布。

参考视频：
中心极限定理（central limit theorem）-统计学
通俗统计学原理入门4 均值抽样分布 中心极限定理

ChatGPT：中心极限定理指出，当从总体（Population）中进行独立随机抽样，并且样本容量（sample size）足够大时，样本均值的分布将接近正态分布，无论总体分布是什么形状。
关于100个平均值的平均值（也称为抽样分布的平均值），它通常被称为抽样均值的均值（mean of sample means）。根据中心极限定理，该抽样均值的均值将趋近于总体的均值。换句话说，通过对多个样本进行抽样并计算平均值，这些样本平均值的平均值将逐渐接近总体的真实均值。

若每次抽样取50个人求平均值，抽100次，这100个平均值的分布仍然会是正态分布，而且mean of sample means还是一样。不同的是，这次的标准差（standard deviation）更小（数据更集中，正态分布的尖更尖）

假设总体工资的标准差为\(\sigma\)，抽样检测得到的“平均值的分布的标准差”（即标准误——standard error，见：通俗理解标准差和标准误)为\(\sigma'\)，每次抽样的样本容量为n，则有如下关系：

\[\frac{\sigma}{\sqrt{n}} = \sigma' \]

第一次抽样每次20人，则$ \sigma_1 = \frac{\sigma}{\sqrt{20}} $
第二次抽样每次100人，则$ \sigma_2 = \frac{\sigma}{\sqrt{100}} $
一般来说， $ \sigma $ 不知道，也可以用某次样本的标准差 \(s\)来代替。
（个人理解：但这样的话应该就不满足正态分布了，而是满足t分布，但只要样本容量足够大，就可以近似认为是正态分布。参考视频：十分钟理解t分布和t检验）
[23-11-10/15:58] 突然发现这篇居然有24的阅读量，还是打个补丁，上面的个人理解是有误的，样本均值的抽样分布符不符合正态是要看样本容量大小的，详细的可以看下文（6.2 总体均值的检验），或者直接参考贾俊平的《统计学》一书

假设检验

已知新冠患者自愈的平均时间是72小时，标准差是8小时。
现在对100人进行试验，发现平均的痊愈时间是69.6小时。
你认为该药物是否有效？
H0：药物无效
H1：药物有效

假设：H0成立
what to do：在H0成立的前提下，计算“出现69.6小时的概率”

因为抽样得到的数据满足正态分布（根据中心极限定理），而且样本平均值也应该是72小时。此时只需看看69.6小时在该正态分布中出现的概率即可。

根据计算得到：在H0成立时，出现69.6小时或更极端的概率是0.3%（P值），我有99.7%的信心拒绝H0。

研究人员的预设：99.5%的信心拒绝H0（置信水平confidence level），（0.5%显著性水平significance level）

参考视频
关于假设检验的一切 - 统计学

t分布

参考视频：
通俗统计学原理入门6 关键一步 从均值抽样分布到t分布
t分布（t distribution）- 统计学

t分布根据我的个人理解，就是中心极限定理，但是把总体的标准差换成样本的标准差。

t检验

t检验我的理解就是一个标准化的过程。不然每次都要去根据标准误算概率，会很蛋疼。

方差分析（ANOVA）

参考视频：
单因素方差分析(上）/ANOVA/什么是方差分析、方差分析的思路

方差分析检验什么——n个分类，它们的某一特征值的平均值，是否有显著差别。

n个分类，比如说B站视频可以分为生活区、鬼畜区、搞笑区……某一特征值，比如说视频播放量。所以方差分析可以分析每个分区的播放量平均值是否相等。

其实就是比较n个样本平均值是否相等，n = 2时直接用t检验即可；n>2的话就用方差分析了。

检验方式：假设检验

\[H_0:\mu_1=\mu_2=\dots =\mu_n \]

\[H_1:\mu_1,\mu_2,\dots ,\mu_n 不全相等 \]

方差分析的思路

数据整体波动 = 组内波动 + 组间波动

数据整体波动（sum of squares total，SST）：B站所有视频播放量的离散程度；

组内波动（sum of squares within，SSW）：搞笑区视频播放量的离散程度；

组间波动（sum of squares between，SSB）：搞笑区、鬼畜区、生活区……播放量的平均值的离散程度；

SSW越大，SSB越小，各组均值相等的可能性越大；
SSW越小，SSB越大，各组均值相等的可能性越小；

方差分析中的计算

参考视频：
单因素方差分析(下)/ANOVA/方差分析的计算、使用Excel和Stata进行方差分析

第5章 参数估计

5.1 参数估计的原理

点估计与区间估计

评量估计量的标准

1. 无偏性（unbiasedness）：
   估计量的抽样分布的期望值等于被估计的总体参数，即$E(\hat{\theta}) = \theta $ （假设\(\theta\)为被估计的总体参数）

由统计量的抽样分布可知，\(E(\bar{x})=\mu, E(p)=\pi, E(s^2)=\sigma^2\)，因此\(\bar{x}、p、s^2\)分别是总体均值\(\mu\)、总体比例\(\pi\)、总体方差\(\sigma^2\)的无偏估计量。

用python进行无偏性模拟：

import numpy as np

# 从均值为50，方差为100的正态总体中，随机抽取10000个样本量为10的样本
# 分别计算出10000个样本均值的均值、样本中位数的均值、样本方差的均值
x, m, v = [], [], []
n = 10

for i in range(1000):
    # loc表示均值，scale表示标准差，size表示生成随机数数量
    d = np.random.normal(loc=50, scale=10, size=n)
    x.append(np.mean(d))
    m.append(np.median(d))
    # np.var计算方差，ddof表示自由度，默认为0，表示总体方差；ddof=1时，表示样本方差；
    v.append(np.var(d, ddof=1))

print(f"样本均值的均值:{np.mean(x):.4f}\n样本中位数的均值:{np.mean(m):.4f}\n样本方差的均值:{np.mean(v):.4f}")

结果为：

样本均值的均值:49.9230
样本中位数的均值:49.9335
样本方差的均值:99.0639

2. 有效性（efficiency）：
   同一总体参数的多个无偏估计量，有更小的方差的估计量更有效。

用python进行有效性模拟：

import numpy as np
# 从均值为0，方差为1的正态总体中，随机抽取10000个样本量为10的样本
# 分别计算出10000个样本均值的方差、样本中位数的方差

x, m = [], []
n = 10
for i in range(10000):
    d = np.random.normal(size=n)
    x.append(np.mean(d))
    m.append(np.median(d))

print(f"样本均值的方差:{np.var(x, ddof=1):.4f}\n样本中位数的方差:{np.var(m, ddof=1):.4f}")

结果为：

样本均值的方差:0.0986
样本中位数的方差:0.1366

绘制样本均值和样本中位数分布的直方图

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False

def plot_dis(df, ax, xlabel):
    df.plot(bins=20, kind='hist', density=True, ax=ax, legend=False)
    df.plot(kind='density', linewidth=2, ax=ax, legend=False)
    ax.set_ylim(0, 1.3)
    ax.set_xlim(-1.5, 1.5)
    ax.set_xlabel(xlabel)
    ax.set_title(xlabel + '的分布')

plt.subplots(1, 2, figsize=(12, 5))
ax1 = plt.subplot(121)
plot_dis(pd.DataFrame(x), ax1, '样本均值')
ax2 = plt.subplot(122)
plot_dis(pd.DataFrame(m), ax2, '样本中位数')

plt.show()

3. 一致性（consistency）：
   随着样本量的增大，统计量收敛于所估计的总体的参数。

样本均值即为总体均值的一个一致估计量。（其实就是随着样本容量增大，样本均值的标准误越来越小——趋向于0——）
（因为\(\sigma_{\bar{x}} = \sigma/\sqrt{n}\)）

用python模拟一致性

import pandas as pd; import numpy as np
# 设置随机种子以再现结果
np.random.seed(2020)
# 生成了一个包含 1000 个随机数的数组 N
N = np.random.normal(loc=50, scale=10, size=1000)
mu = np.mean(N)

# np.random.choice() 函数用于从给定的数组中进行随机抽样。
# 它接受三个参数：抽样的数据源（这里是数组 N）、抽样的数量（这里是 10）、是否允许重复抽样（这里设置为 replace=False 表示不允许重复抽样）。
xbar10 = np.mean(np.random.choice(N, 10, replace=False))
xbar100 = np.mean(np.random.choice(N, 100, replace=False))
xbar500 = np.mean(np.random.choice(N, 500, replace=False))
xbar900 = np.mean(np.random.choice(N, 900, replace=False))

# .T表示转置
aver = pd.DataFrame([mu, xbar10, xbar100, xbar500, xbar900],
                   ['总体均值', 'xbar10', 'xbar100', 'xbar500', 'xbar900']).T

print(aver)

# 样本均值与总体均值mu的差值d
d = pd.DataFrame([{'d10':(xbar10-mu), 'd100': (xbar100-mu), 'd500': (xbar500-mu), 'd900': (xbar900-mu)}])
print(d)

结果为：

   总体均值    xbar10    xbar100    xbar500    xbar900
0  49.665586  48.45996  50.165167  49.568964  49.636309
     d10      d100      d500      d900
0 -1.205626  0.499581 -0.096622 -0.029277

5.2 总体均值的区间估计

研究一个总体，推断总体均值\(\mu\)的统计量就是样本均值\(\bar{x}\)。
研究两个总体时，所关心的参数主要是两个总体均值之差\((\mu_1 - \mu_2)\)，用于推断的统计量则是两个样本的均值之差$ (\bar{x_1} - \bar{x_2})$。

5.2.1一个总体均值的估计

要分大样本\(n \ge 30\) 和小样本$ n < 30 $ 两种情况讨论。
不管哪种情况，总体均值的置信区间都是由样本均值加减估计误差得到的。
估计误差由（1）点估计量的标准误；（2）估计时所要求的置信水平为\((1-\alpha)\)时统计量分布两侧面积各为\(\alpha/2\)时的分位数值。

总体均值在\((1-\alpha)\)置信水平下的置信区间可一般性地表达为：

\[\bar{x} \pm (分位数 \times \bar{x}的标准误) \]

1. 大样本的估计

由中心极限定理，样本均值\(\bar{x}\) 经过标准化后服从标准正态分布，即

\[z = \frac{\bar{x}-\mu}{\sigma/\sqrt{n}} \sim N(0,1) \]

当总体标准差\(\sigma\)未知时，可以用样本标准差\(s\)代替。

当总体标准差\(\sigma\)已知时，总体均值\(\mu\)在\((1-\alpha)\)置信水平下的置信区间为：

\[\bar{x} \pm z_{\alpha/2}\frac{\sigma}{\sqrt{n}} \]

当总体标准差\(\sigma\)未知时，总体均值\(\mu\)在\((1-\alpha)\)置信水平下的置信区间为：

\[\bar{x} \pm z_{\alpha/2}\frac{s}{\sqrt{n}} \]

式中$ \bar{x} - z_{\alpha/2}\frac{\sigma}{\sqrt{n}}$为置信下限；

$ \bar{x} + z_{\alpha/2}\frac{\sigma}{\sqrt{n}}$为置信上限；

\(\alpha\)为事先确定的一个概率值，它是总体均值不包含在置信区间内的概率；

\((1-\alpha)\)为置信水平；

\(z_{\alpha/2}\)为标准正态分布两侧面积各位\(\alpha/2\)时的\(z\)值；

\(z_{\alpha/2}\frac{\sigma}{\sqrt{n}}\)为估计误差。

用python模拟：

import pandas as pd
import numpy as np
import scipy.stats as st

# 40辆车每百公里的耗油量数据
example5_1 = pd.read_csv('D:/DataDisk/307923JiaJunPing_data_analysis/example/chap05/example5_1.csv', encoding='gbk')

# scipy.stats.norm.interval 函数的作用是根据给定的置信水平和正态分布的均值与标准差，计算出对应的置信区间。
# 置信区间表示了对总体参数（例如均值）的估计范围，通常以置信水平的形式给出（例如 95% 置信水平）。
# scipy.stats.sem(arr，axis = 0，ddof = 0)函数用于计算输入数据平均值的标准误差。（应该就是标准差吧）
int = st.norm.interval(0.90, loc=np.mean(example5_1), scale=st.sem(example5_1))
# 以数组形式返回置信区间
int_4 = np.round(int, 4)

print(f"置信区间为：{int_4}")
print(f"耗油量平均值为：{np.mean(example5_1)}")
print(f"耗油量标准差为：{st.sem(example5_1)}")

结果为：

置信区间为：[[7.8359]
 [8.0991]]
耗油量平均值为：7.967500000000001
耗油量标准差为：[0.08001502]

2. 小样本的估计

小样本情形下，需要假设总体服从正态分布。

在估计之前首先应对总体的正态性进行检验，可以用直方图、茎叶图、P-P图或Q-Q图进行初步评估，也可以使用Shapiro检验和K-S检验等，具体内容见第6章。

• 若正态总体的标准差已知，样本均值服从正态分布：

（这种情形前面已经讨论过了）

• 若正态总体的标准差未知，用样本标准差代替总体标准差，样本均值服从t分布：

在\((1-\alpha)\)置信水平下，总体均值的置信区间为：

\[\bar{x} \pm t_{\alpha/2}\frac{s}{\sqrt{n}} \]

import pandas as pd
import numpy as np
import scipy.stats as st

# 25袋食品重量的数据
example5_2 = pd.read_csv('D:/DataDisk/307923JiaJunPing_data_analysis/example/chap05/example5_2.csv', encoding='gbk')

# 第二个参数是自由度
int = st.t.interval(0.95, len(example5_2)-1, loc=np.mean(example5_2), scale=st.sem(example5_2))

int_4 = np.round(int, 4)

print(f"食品平均质量的置信区间为：{int_4}")
print(f"食品重量平均值为：{np.mean(example5_2)}")
print(f"食品重量标准差为：{st.sem(example5_2)}")

结果为：

食品平均质量的置信区间为：[[101.3748]
 [109.3452]]
食品重量平均值为：105.36000000000001
食品重量标准差为：[1.93089789]

5.2.2 两个总体均值差的估计

第6章 假设检验

反转了，贾俊平主编的数据分析基础——基于python的实现一书中，又说大样本的情况下，不知道总体均值也可以用样本均值代替，因为样本是大样本。

6.2 总体均值的检验

6.2.1 一个总体均值的检验

1. 大样本的检验

在大样本（n>30）的情形下，样本均值的抽样分布服从正态分布，其标准误为$ \sigma/\sqrt{n} $ ；若总体均值为\(\mu_0\)，总体方差\(\sigma^2\)已知时，可以进行z检验：

\[z = \frac{\bar{x}-\mu_0}{\sigma/\sqrt{n}} \]

当总体方差\(\sigma^2\)未知时，可以用样本方差\(s^2\)来代替：

\[z = \frac{\bar{x}-\mu_0}{s/\sqrt{n}} \]

【例6-3】城市过去每立方米空气中PM2.5的均值是\(81\mu g/m^3\)。根据最近的40次检验（example6_3.csv），能否认为城市的PM2.5均值显著低于\(81\mu g/m^3\)？

• 假设（这是单侧检验）：
  \(H_0: \mu \ge 81\)
  \(H_1: \mu < 81\)

• 检验（用python）

import pandas as pd
from statsmodels.stats.weightstats import ztest

example6_3 = pd.read_csv('path/to/data/example6_3.csv', encoding="gbk")
z, p_value = ztest(x1=example6_3['PM2.5值'], value=81, alternative='smaller')

print(f"样本均值 = {example6_3['PM2.5值'].mean():.2f}, z统计量值 = {z:.4f}, p值 = {p_value:.4f}")

输出结果为：

样本均值 = 79.55, z统计量值 = -1.1856, p值 = 0.1179

由于\(P>0.05\)，相信原假设。

2. 小样本的检验

在小样本（n<30）的情形下，首先假定总体服从正态分布。

所以如果总体不服从正态分布呢？书里没说。留个坑。

[23-10-30] 看了贾俊平的《统计学》，里面说的是，小样本的情形下，如果总体方差已知，那么样本均值服从正态分布，用z检验；如果总体方差未知，用样本方差代替，则用t检验。

• 当总体方差\(\sigma^2\)已知时，样本均值服从正态分布，可以用z检验。

• 当总体方差未知时，用样本方差\(s^2\)代替\(\sigma^2\)，此时样本均值$ \sim t(n-1) $，可以用t检验。

\[t = \frac{\bar{x}-\mu_0}{s/\sqrt{n}} \]

效应量：表示样本均值与假设的总体均值相差多少个标准差。

在单样本t检验中通常使用Cohen的\(d\) 统计量来度量：

\[d = \frac{|样本均值-假设的总体均值|}{样本标准差} = \frac{|\bar{x}-\mu_0|}{s} \]

根据Cohen（1988）提出的标准：

小效应量	中效应量	大效应量
0.20	0.50	0.80

$ d > 0.80 $时，为大的效应量，表示相差0.8个标准差。

【例6-4】砖头的厚度要求为5cm。现抽取20块砖进行检验。
假定砖的厚度服从正态分布，在0.05的显著性水平下，检验砖的厚度是否合格。

• 假设：
  \(H_0: \mu = 5\)
  \(H_1: \mu \ne 5\)

• （用python）检验：

import pandas as pd
from scipy.stats import ttest_1samp

example6_4 = pd.read_csv('path/to/data/example6_4.csv', encoding="gbk")

t, p_value = ttest_1samp(a=example6_4['厚度'], popmean=5)  # popmean为假设的总体均值
print(f"样本均值 = {example6_4['厚度'].mean():.2f}, t统计量值 = {t:.4f}, p值 = {p_value:.4g}")

输出结果为：

样本均值 = 4.80, t统计量值 = -5.6273, p值 = 1.998e-05

\(P<0.05\) 放弃原假设，有证据显示砖头厚度不合格。

• （用python）计算效应量：

example6_4 = pd.read_csv('path/to/data/example6_4.csv', encoding="gbk")
mu = 5
xbar = example6_4['厚度'].mean()
sd = example6_4['厚度'].std()
d = abs(xbar-mu)/sd
print(f"效应量d={d:.4f}")

输出结果为：

效应量d=1.2583

表明样本的砖头平均厚度与标准厚度相差1.2583个标准差，根据Cohen准则，属于比较大的效应量。

6.2.2 两个总体均值差的检验

根据获得样本的方式不同，分为独立样本和配对样本两种情形。

1. 独立大样本检验

在大样本情形下，两个样本均值之差 $ \bar{x}_1 - \bar{x}_2 $ 的抽样分布近似服从正态分布。

故可用z检验。

• 若两个总体的方差\(\sigma_1\)、\(\sigma_2\)已知，则：

\[z = \frac{(\bar{x}_1 - \bar{x}_2)-(\mu_1-\mu_2)}{\sqrt{\frac{\sigma_1^2}{n_1} - \frac{\sigma_2^2}{n_2}}} \]

• 若两个总体的方差\(\sigma_1\)、\(\sigma_2\)未知，用样本方差$ s_1^2 \(、\)s_2^2$ 代替，则：

\[z = \frac{(\bar{x}_1 - \bar{x}_2)-(\mu_1-\mu_2)}{\sqrt{\frac{s_1^2}{n_1} - \frac{s_2^2}{n_2}}} \]

用python模拟：

import pandas as pd
from statsmodels.stats.weightstats import ztest

# 分析男女学生上网时间是否有差异
# 原假设：男女学生上网时间没有差异，即μ1 - μ2 = 0
# 文件为男生上网时间（第一列）和女生上网时间（第二列）
example6_5 = pd.read_csv('path\to\data\chap06\example6_5.csv', encoding="gbk")

# alternative='two-sided'表示双侧检验
z, p_value = ztest(x1=example6_5['男生上网时间'], x2=example6_5['女生上网时间'], alternative='two-sided')
print(f"男生平均上网时间为 = {example6_5['男生上网时间'].mean():.4f}\
\n女生平均上网时间为 = {example6_5['女生上网时间'].mean():.4f}\
\nz统计值 = {z:.4f} \
\np值 = {p_value:.4f}")

输出结果为：

男生平均上网时间为 = 3.0583
女生平均上网时间为 = 2.8306
z统计值 = 1.1188
p值 = 0.2632

表明无证据显示男女学生上网的平均时间有显著差异。

2. 独立小样本的检验

首先假设两个总体均服从正态分布。

分三种情形：

1. 两个正态总体方差已知——无论样本大小，两个样本均值之差的抽样分布都服从正态分布。（上文已讨论过检验方法）

2. 两个正态总体的方差未知，但已知\(\sigma_1^2 = \sigma_2^2\)。

因为两个总体的方差相等，故可以把两个样本的方差合并起来，以估计总体方差。

这个被合并起来的估计量叫做总体方差的合并估计量（pooled variance/pooled sample variance），记为\(s_p^2\)。

\[\frac{s_p^2 = (n_1 - 1)s_1^2+(n_2-1)s_2^2}{n_1+n_2 -2} \]

（计算出来的值将介于\(s_1^2\)和\(s_2^2\)之间，但会偏向于样本容量较大的那个样本方差）

这时，两个样本均值之差经标准化后服从自由度为\((n_1+n_2-2)\) 的t分布，故采用t检验：

\[t = \frac{(\bar{x}_1 - \bar{x}_2)-(\mu_1-\mu_2)}{s_p\sqrt{\frac{1}{n_1} - \frac{1}{n_2}}} \]

参考视频：Pooled-Variance t Tests and Confidence Intervals: Introduction

3. 两个正态总体的方差未知且不相等时，两个样本均值之差经标准化后近似服从自由度为\(v\)的t分布。

\[t = \frac{(\bar{x}_1 -\bar{x}_2)-(\mu_1-\mu_2)}{\sqrt{\frac{s_1^2}{n_1} - \frac{s_2^2}{n_2}}} \]

\(t\)的自由度\(v\)为：

\[v = \frac{(\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2})^2}{ \frac{(s_1^2/n_1)^2}{n_1-1} + \frac{(s_2^2/n_2)^2}{n_2-1}} \]

用python模拟：

import pandas as pd
from statsmodels.stats.weightstats import ttest_ind


# 分析两家企业生产的灯泡的平均使用寿命是否有显著差异
# 文件为甲企业（第一列）和乙企业（第二列）数据
example6_6 = pd.read_csv('path\to\data\chap06\example6_6.csv', encoding="gbk")

xbar1 = example6_6['甲企业'].mean()
xbar2 = example6_6['乙企业'].mean()

# （1）假设总体方差相等
# pooled表示总体方差相等，unequal为不等
t, p_value, df = ttest_ind(x1=example6_6['甲企业'], x2=example6_6['乙企业'], 
                           alternative='two-sided',
                           usevar='pooled')
print("==============（1）假设总体方差相等==============")
print(f"甲企业的灯泡平均寿命 {xbar1}\
\n乙企业的灯泡平均寿命 = {xbar2}\
\nt统计值 = {t:.6f} \
\n自由度 = {df} \
\np值 = {p_value:.6f}")

# （2）假设总体方差不等

t, p_value, df = ttest_ind(x1=example6_6['甲企业'], x2=example6_6['乙企业'], 
                           alternative='two-sided',
                           usevar='unequal')
print("==============（2）假设总体方差不等==============")
print(f"甲企业的灯泡平均寿命 {xbar1}\
\n乙企业的灯泡平均寿命 = {xbar2}\
\nt统计值 = {t:.6f} \
\n自由度 = {df:.4f} \
\np值 = {p_value:.6f}")

结果为：

==============（1）假设总体方差相等==============
甲企业的灯泡平均寿命 8487.5
乙企业的灯泡平均寿命 = 8166.0
t统计值 = 3.494270
自由度 = 38.0
p值 = 0.001225
==============（2）假设总体方差不等==============
甲企业的灯泡平均寿命 8487.5
乙企业的灯泡平均寿命 = 8166.0
t统计值 = 3.494270
自由度 = 33.6826
p值 = 0.001353

结果表明两家企业的灯泡平均寿命有显著区别。

差别的程度由效应量给出。

独立样本t检验的效应量的估计通常由Cohen的d统计量给出，计算公式为：

\[d = |t|\sqrt{\frac{n_1+n_2}{n_1n_2}} \]

# 承接前文代码，沿用t值
n1 = len(example6_6['甲企业'])
n2 = len(example6_6['乙企业'])

d = abs(t)*np.sqrt((n1+n2)/(n1*n2))

print(f"效应量d={d:.6f}")

效应量d=1.104985

3. 配对样本的检验

假设两个总体配对的差值服从正态分布，而且配对差是从差值总体中随机抽取的。

• 对于大样本情形，配对差值服从正态分布，进行z检验。

• 对于小样本情形，配对差值服从t分布，自由度为n-1，进行t检验。

\[t = \frac{\bar{d} - (\mu_1 - \mu_2)}{s_d/\sqrt{n}} \]

\(\bar{d}\) 为配对差值的均值；\(s_d\)为配对差值的标准差。

用python模拟：

import pandas as pd
from scipy.stats import ttest_rel

# 分析消费者对两款饮料的评分是否有差异
# H0：评分无差异，即μ1 - μ2 = 0 
# 文件为旧款饮料（第一列）和新款饮料（第二列）的评分，样本量为10
example6_7 = pd.read_csv('path\to\data\chap06\example6_7.csv', encoding="gbk")

dbar = (example6_7['旧款饮料'] - example6_7['新款饮料']).mean()

t, p_value = ttest_rel(a = example6_7['旧款饮料'], b=example6_7['新款饮料'])

print(f"配对样本差值的均值 {dbar}\
\nt统计值 = {t:.6f} \
\np值 = {p_value:.6g}")

配对样本差值的均值 -1.3
t统计值 = -2.750848
p值 = 0.0224457

结果表明有显著差异。

计算效应量：

\[d = \frac{|配对样本差值的均值|}{配对差值的标准差} = \frac{\bar{d}}{s_d} = \frac{|t|}{\sqrt{n}} \]

n = example6_7.shape[0]

d = abs(t)*np.sqrt(n)

print(f"效应量d={d:.6f}")

效应量d=8.698945

第8章 方差分析

8.1 方差分析的原理

原理还是看前面《统计基础》的笔记吧，这里列一些专用名词。

例子：不同视频分区（如：生活区、舞蹈区、鬼畜区、学习区、游戏区）的播放量是否有显著区别。

名词	英文	解释	对应例子
类别变量	categorical variable	只能取有限可能值的变量	性别
因子	factor	方差分析中的类别变量	B站视频的分区
处理/水平	treatment/level	因子的取值	生活区、舞蹈区……
试验单元	experiment unit	（这不知咋解释）	每个B站视频

---

名词	英文	解释
总平方和	sum of squares for total（SST）	反映全部数据总误差大小
处理平方和/组间平方和	treatment sum of squares（SSA，这里把因子记为A）/between-group sum of squares	其实就是SSB
误差平方和/组内平方和	sum of squares of error（SSE）/within-group sum of squares	其实就是SSW

8.2 单因子方差分析（one-way analysis of variance)

数学模型

设因子A有I个处理，则

\[y_{ij} = \mu_i + \epsilon_{ij} \]

\(y_{ij}\)为第\(i\)个处理中的第\(j\)个观测值；

\(\mu_i\)为第\(i\)个处理的平均值；

\(\epsilon_{ij}\)为第\(i\)个处理的第\(j\)个观测值的随机误差。

设全部观测数据的总均值为\(\mu\)，第i个处理的处理效应为\(\alpha_i = \mu_i-\mu\)，则前面的式子可写为：

\[y_{ij} = \mu + \alpha_i + \epsilon_{ij} \]

• 假设：
  $H_0: \alpha_1 = \alpha_2 = \dots =\alpha_I = 0 \(（不同处理没有显著效应） \)H_1: \(至少有一个\)\alpha_i \ne 0 $（不同处理有显著效应）

等价的假设形式为：
$H_0: \mu_1 = \mu_2 = \dots =\mu_I $
$H_1: \mu_i $们不全相等

单因子方差分析表：

	平方和SS	自由度df	均方MS	检验统计量
处理效应	$$SSA=\sum_{i=1}^I n_i(\bar{y_i}-\bar{\bar{y}})^2$$	$$I-1$$	$$MSA=\frac{SSA}{I-1}$$	$$\frac{MSA}{MSE}$$
误差	$$SSE = \sum_{i=1}^I \sum_{j=1}^{n_i} (y_{ij}-\bar{y_i})^2$$	$$ n-I$$	$$MSE=\frac{SSE}{n-I}$$
总效应	$$ SST=\sum_{i=1}^I \sum_{j=1}^{n_i} (y_{ij}-\bar{\bar{y}})^2 $$	$$n-1$$

n为观测值个数——B站视频的个数；

\(n_i\)为第i个处理的样本容量——舞蹈区的视频总数；

\(\bar{y_i}\)为第i个处理的样本均值——舞蹈区视频的播放量的平均值；

\(\bar{\bar{y}}\) 是所有视频的播放量的平均值；

求出F值后，就可以进行F检验——计算P值。

用python模拟

import pandas as pd
from statsmodels.formula.api import ols
from statsmodels.stats.anova import anova_lm

# 文件为3个不同季度的CPU性能测试得分表
example8_1 = pd.read_csv('path/to/data/example8_1.csv', encoding='gbk')

# 这行的作用就是把数据列表的列名给换了
# 第一列列名换成season，第二列……
example8_1.columns = ['season', 'value']

# 这段代码是在使用OLS（Ordinary Least Squares，普通最小二乘法）进行线性回归分析。
# 首先，'value~C(season)' 是一个公式字符串，它指定了线性回归模型的形式。
# 具体来说，它表示依据 'value' 列对 'season' 列进行回归分析。
# C(season) 表示将 'season' 列视为分类变量（categorical variable），而不是连续变量。
model = ols('value~C(season)', example8_1[['season', 'value']]).fit()

# 对线性回归模型进行方差分析
a = anova_lm(model)
print(a)

输出结果为

             df        sum_sq       mean_sq          F        PR(>F)
C(season)   2.0  4.590439e+12  2.295219e+12  49.708864  4.007688e-15
Residual   87.0  4.017072e+12  4.617324e+10        NaN           NaN

检验的P值（PR(>F)）= 4.007688e-15<0.05，舍弃假设H0，意味着评测时间对CPU性能有影响

至于在代码中为什么要先做线性回归分析，再进行方差分析，还没搞清楚。

搞清楚了，好像是因为，方差分析其实就是一种特殊的广义线性回归分析（参考：Introduction to general linear models

第9章 一元线性回归

9.1 确定变量间的关系

1. Pearson相关系数

相关系数（correlation coefficient）：度量两个变量之间线性关系强度的统计量。记为\(r\)，则

\[r = \frac{\sum(x_i- \bar x )(y_i-\bar y)}{\sqrt{\sum(x_i- \bar x )^2\cdot \sum(y_i-\bar y)^2}} \]

\(r = 1\)表示完全正相关，\(r = -1\)表示完全负相关，$ r = 0 $表示不存在线性关系（但不代表不存在其他关系）

2. 相关系数的检验

总体相关系数\(\rho\)，通常是未知的，需用样本相关系数\(r\)进行估计。

验证相关系数\(r\)的可靠性：t检验方法（R.A.Fisher提出）

1. 提出假设：

\(H_0: \rho = 0（两个标量的线性关系不显著）\)
\(H_1: \rho \ne 0（两个标量的线性关系不显著）\)

2. 计算t值：

\[t = \frac{r\sqrt{n-2}}{\sqrt{1-r^2}} \sim t(n-2) \]

n是什么？留个坑

3. 决策：
   求出\(P\)值，若$P < \alpha \(，拒绝\)H_0$，表示总体的两个变量之间线性关系显著。

用python检验：

import pandas as pd
from scipy.stats import pearsonr

example9_1 = pd.read_csv('path/to/data/example9_1.csv', encoding='gbk')
corr, p_value = pearsonr(example9_1['居民家庭消费支出'], example9_1['居民家庭总收入'])
# pearsonr(x,y)
print(corr, p_value)

9.2 模型估计与检验

回归模型与回归方程

回归模型：

\[y = \beta_0 + \beta_1x+\epsilon \]

回归方程：

\[\hat{y} = \hat{\beta_0} + \hat{\beta_1}x+\epsilon \]

？？？一模一样？，没有帽子的回归模型表示你觉得总体就是这么个情况，其中的\(\beta\)是真实确定有这么个值，只是你不知道大小；有帽子的回归方程表示这是你通过样本信息推测出来的公式，带帽子的值表示都是你根据样本信息推测出来的值。

\(\epsilon\)表示误差项，即不能为模型所解释的偏差。

参数估计

最小二乘估计（least squares estimation）：

\[\hat{\beta_1} = \frac{\sum{(x_i-\bar{x})(y_i-\bar{y})}}{\sum{(x_i-\bar{x})^2}} \]

\[\hat{\beta_0} = \bar{y} - \hat{\beta_1}\bar{x} \]

用python进行参数估计：

import pandas as pd
from statsmodels.formula.api import ols

example9_1 = pd.read_csv('path/to/data/example9_1.csv', encoding='gbk')
model = ols("居民家庭消费支出~居民家庭总收入", data=example9_1).fit()
# old("y~x", data=data).fit()

print(model.summary())

得到输出结果：

                            OLS Regression Results                            
==============================================================================
Dep. Variable:       居民家庭消费支出   R-squared:                       0.980
Model:                            OLS   Adj. R-squared:                  0.980
Method:                 Least Squares   F-statistic:                     2890.
Date:                Sat, 21 Oct 2023   Prob (F-statistic):           3.52e-51
Time:                        18:04:32   Log-Likelihood:                -514.68
No. Observations:                  60   AIC:                             1033.
Df Residuals:                      58   BIC:                             1038.
Df Model:                           1                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept    1.894e+04    643.529     29.436      0.000    1.77e+04    2.02e+04
居民家庭总收入  0.4725      0.009     53.758      0.000       0.455       0.490
==============================================================================
Omnibus:                        6.328   Durbin-Watson:                   2.139
Prob(Omnibus):                  0.042   Jarque-Bera (JB):                2.477
Skew:                           0.102   Prob(JB):                        0.290
Kurtosis:                       2.026   Cond. No.                     2.79e+05
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 2.79e+05. This might indicate that there are
strong multicollinearity or other numerical problems.

估计方程为：

\[\hat{y} = 18942.9031 + 0.4725 \times 居民家庭总收入 \]

模型的拟合优度（goodness of fit）

拟合优度：回归直线与观测点的接近程度
拟合优度的度量方式：决定系数（coefficient of determination）

1. 决定系数

每个观测值的误差都可以分解为两个部分：$ y_i - \bar{y} = (y_i - \hat{y_i})+(\hat{y_i} - \bar{y})（画图特好理解）$

两边平方，并对所有n个点求和，有：

\[\sum{(y_i-\bar{y})^2} = \sum{(y_i-\hat{y_i})^2} + \sum{(\hat{y_i}-\bar{y})^2} + 2\sum{(y_i-\hat{y_i})(\hat{y_i}-\bar{y})} \]

可以证明：\(\sum{(y_i-\hat{y_i})(\hat{y_i}-\bar{y})} = 0\)，因此有

\[\sum{(y_i-\bar{y})^2} = \sum{(y_i-\hat{y_i})^2} + \sum{(\hat{y_i}-\bar{y})^2} \]

等式左边为总平方和SST（total sum of squares）；

右边第一项为残差平方和SSE（residual sum of squares），是不能被回归方程解释的误差部分；

右边第二项为回归平方和SSR（regression sum of squares），此误差是由x的变化引起的，是可以被回归方程解释的部分；

SST = SSE + SSR

决定系数：

\[R^2 = \frac{SSR}{SST} = \frac{\sum{(\hat{y_i}-\bar{y})^2}}{\sum{(y_i-\bar{y})^2}} \]

若所有观测点都落在回归直线上，\(SSR=0\)，\(R^2=1\)，模型是完全拟合的。

\(R^2 \rightarrow 1\) 拟合得越好
\(R^2 \rightarrow 0\) 拟合得越差

在前面的代码示例中，R-squares=0.938=93.8%表示在居民家庭消费支出取值的总误差中，有93.8%可以由居民家庭消费支出与家庭总收入之间的线性关系来解释。

2. 残差的标准误（residual standard error）

残差的标准误（\(s_e\)）是SSE的均方根：

\[s_e = \sqrt{\frac{\sum{(y_i-\hat{y_i})^2}}{n-k-1}} = \sqrt{\frac{SSE}{n-k-1}} = \sqrt{MSE} \]

其中，k为自变量的个数，在一元线性回归中，\(k = 1, n-k-1 = n-2\)

\(s_e\)反映了\(y_i\)在回归直线周围的离散程度。

模型的显著性检验

什么？前面不是已经做过Pearson相关系数检验过线性关系了吗？
问了ChatGPT，大概意思是说，Pearson相关系数检验告诉我们变量之间有线性关系，而接下来的F检验和t检验告诉我们这个模型拟合的好坏。

[Log-2023-11-1]：我现在的理解是，线性关系检验是检验y与{xi}的线性关系显不显著，而回归系数的检验是检验y与具体每个xi的线性关系显不显著。

1. 线性关系检验（F检验）

F检验 检验因变量y和自变量x之间的线性关系是否显著（或者说，它们之间能否用\(y=\beta_0+\beta_1x+\epsilon\)来表示。）

\[F = \frac{SSR/k}{SSE/n-k-1} = \frac{MSR}{MSE} \sim F(k, n-k-1) \]

MSR：回归均方
MSE：残差均方

检验步骤：

• 提出假设：
  $ H_0:\beta_1=0（两个变量之间的线性关系不显著）$
  $ H_1:\beta_1 \ne 0（两个变量之间的线性关系显著）$

• 计算F

• 求F的P值，若 $ P < \alpha $，则拒绝 $ H_0 $

若H0成立，\(MSR/MSE \rightarrow 1\)；
若H0不成立，\(MSR/MSE \rightarrow \inf\)；
MSR/MSE较大时拒绝H0，继而断定x和y之间线性关系显著；

2. 回归系数的检验和判断（t检验）

t检验 检验自变量对因变量的影响是否显著。

在一元线性回归中，由于只有一个自变量，回归系数检验与线性关系检验是等价的
（在多元线性回归中不再等价）

检验步骤：

• 假设：
  \(H_0: \beta_1 = 0 （自变量对因变量的影响不显著）\)
  \(H_1: \beta_1 \ne 0 （自变量对因变量的影响显著）\)

检验统计量的构造是以回归系数\(\beta_1\)的抽样分布为基础的（这句看不懂，所以照抄下来）

统计证明，\(\hat{\beta_1}\)服从正态分布，期望值为\(E(\hat{\beta_1})=\beta_1\)，标准误的估计量为：

\[s_{\hat{\beta_1}} = \frac{s_e}{\sqrt{\sum{x_i^2 - \frac{1}{n}(x_i)^2}}} \]

（\(s_e\)为残差标准误，在前文中有公式）

将回归系数标准化，即可得到检验\(\beta_1\)的统计量t。

在H0成立的条件下，\(\hat{\beta_1}-\beta_1 = \hat{\beta_1}\)

故：

\[t = \frac{\hat{\beta_1}}{s_{\hat{\beta_1}}} \sim t(n-2) \]

• 计算t值

• 求t的P值，若 $ P < \alpha $ ，则拒绝 $ H_0 $ ，表示x对y的影响显著。

第10章 多元线性回归

10.1 多元线性回归模型及其参数估计

回归方程与回归模型

参数的最小二乘估计

使残差平方和最小，即：

\[Q = \sum{(y_i-\hat{y_i})^2} = \sum{(y_i- \hat{\beta_0} - \hat{\beta_1}x_1 - \dots - \hat{\beta_k}x_k)^2} = min \]

得到求解$ \hat{\beta_0}, \hat{\beta_1}, \dots, \hat{\beta_k} $的标准方程组为：

\[\left. \frac{\partial Q}{\partial \beta_0} \right|_{\beta_0 = \hat{\beta_0}} = 0 \]

\[\left. \frac{\partial Q}{\partial \beta_i} \right|_{\beta_i = \hat{\beta_i}} = 0, i = 1, 2, \dots, k \]

用python估计：

# 多元线性回归
import pandas as pd
from statsmodels.formula.api import ols

example10_1 = pd.read_csv('path/to/data/example10_1.csv', encoding='gbk')
model = ols("y~x1+x2+x3+x4", data=example10_1).fit()
# y为居民家庭消费支出；x1为居民家庭总收入；
# x2为受访者年龄；x3为受访者全年总收入；x4为家庭常住人口；
print(model.summary())

得到输出结果：

                            OLS Regression Results                            
==============================================================================
Dep. Variable:                      y   R-squared:                       0.985
Model:                            OLS   Adj. R-squared:                  0.984
Method:                 Least Squares   F-statistic:                     888.2
Date:                Sat, 21 Oct 2023   Prob (F-statistic):           3.05e-49
Time:                        18:52:46   Log-Likelihood:                -507.03
No. Observations:                  60   AIC:                             1024.
Df Residuals:                      55   BIC:                             1035.
Df Model:                           4                                         
Covariance Type:            nonrobust                                         
==============================================================================
                 coef    std err          t      P>|t|      [0.025      0.975]
------------------------------------------------------------------------------
Intercept   2.438e+04   1534.983     15.882      0.000    2.13e+04    2.75e+04
x1             0.4532      0.013     36.149      0.000       0.428       0.478
x2           -56.6435     18.228     -3.107      0.003     -93.173     -20.114
x3            -0.0097      0.010     -0.943      0.350      -0.030       0.011
x4          -343.1776    150.159     -2.285      0.026    -644.103     -42.252
==============================================================================
Omnibus:                        1.715   Durbin-Watson:                   2.374
Prob(Omnibus):                  0.424   Jarque-Bera (JB):                1.275
Skew:                           0.116   Prob(JB):                        0.529
Kurtosis:                       2.325   Cond. No.                     8.98e+05
==============================================================================

Notes:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
[2] The condition number is large, 8.98e+05. This might indicate that there are
strong multicollinearity or other numerical problems.

估计方程为：

\[\hat{y} = 24378.2034 + 0.4532x_1 - 56.6435x_2 - 0.0097 x_3 - 343.1776x_4 \]

第11章 时间序列分析

11.1 时间序列的成分和预测方法

记号：

\(t\) 表示所观察的时间
\(\hat{Y_t}(t=1,2, \dots , n)\) 表示在时间 \(t\) 上的观测值。

时间序列的成分

1. 趋势（trend）；
2. 季节变动（seasonal fluctuation）：以年为周期的固定变动；
3. 循环波动（cyclical fluctuation）：非固定长度的周期性波动；
4. 不规则波动（irregular variations）；

预测方式

预测方法	适合的数据模式	对数据的要求	预测期
简单指数平滑	随机波动	5个以上	短期
Holt指数平滑	线性趋势	5个以上	短期至中期
一元线性回归	线性趋势	10个以上	短期至中期
指数模型	非线性趋势	10个以上	短期至中期
多项式函数	非线性趋势	10个以上	短期至中期
Winters指数平滑	趋势和季节成分	至少有4个周期的季度或月份数据	短期至中期
分解预测	趋势、季节和循环成分	至少有4个周期的季度或月份数据	短、中、长期

评估方法的好坏：

• 平均误差（mean error）
• 平均绝对误差（mean absolute deviation）
• 均方误差（mean square error）
• 平均百分比误差（mean percentage error）
• 平均绝对百分比误差（mean absolute percentage error）

常用：均方误差；
均方误差（MSE）是误差平方和的平均数

\[MSE = \frac{\sum_{i=1}^{n}(Y_i-F_i)^2}{n} \]

\(Y_i\) 为第i期的实际值；\(F_i\) 为第i期的预测值；$ n $ 为预测误差的个数；