预测（Prediction）

除了探究数据本身的性质，数据科学的另一个重点就是基于已有数据进行预测。最常见的预测就是基于一个变量预测另一个变量的值。

在之前的身高预测例子中，我们通过观察散点图得到变量之间的关联性，并得到了一条回归线。下面我们会在此基础上通过衡量样本与回归线的距离进一步刻画变量之间的关系。

实际上，预测与推断现代统计学/机器学习语义下已经没有明显的界限，但在本课程中，“推断”还是专指统计学上的参数估计与假设检验，而预测则表示对一个具体值的推算。

（线性）相关性（Correlation）

首先我们探究散点图中数据围绕回归线分布的紧密程度，即线性关系。以下面这个hybrid数据集为例（前10行）：

车型	年份	建议售价	加速度	油耗	类别
Prius (1st Gen)	1997	24509.7	7.46	41.26	Compact
Tino	2000	35355	8.2	54.1	Compact
Prius (2nd Gen)	2000	26832.2	7.97	45.23	Compact
Insight	2000	18936.4	9.52	53	Two Seater
Civic (1st Gen)	2001	25833.4	7.04	47.04	Compact
Insight	2001	19036.7	9.52	53	Two Seater
Insight	2002	19137	9.71	53	Two Seater
Alphard	2003	38084.8	8.33	40.46	Minivan
Insight	2003	19137	9.52	53	Two Seater
Civic	2003	14071.9	8.62	41	Compact

为了衡量变量之间的线性关系，我们需要先对变量进行标准化。最常用的标准化就是Z-score标准化，即减去均值后除以标准差，示例：

import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

hybrid = pd.read_csv('hybrid.csv')
suv = hybrid[hybrid['class'] == 'SUV']

# Z-score 函数
def standard_units(any_numbers):
    return (any_numbers - np.mean(any_numbers)) / np.std(any_numbers)

suv['mpg_std'] = standard_units(suv['mpg'])
suv['msrp_std'] = standard_units(suv['msrp'])
suv['acceleration_std'] = standard_units(suv['acceleration'])  

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), dpi=300)

# mpg vs msrp
ax1.scatter(suv['mpg_std'], suv['msrp_std'])
ax1.set_xlim(-3, 3)
ax1.set_ylim(-3, 3)
ax1.set_xlabel('mpg (standard units)')
ax1.set_ylabel('msrp (standard units)')
ax1.grid(True, linestyle='--', alpha=0.7)
ax1.set_title('MPG vs MSRP')

# acceleration vs msrp
ax2.scatter(suv['acceleration_std'], suv['msrp_std'])
ax2.set_xlim(-3, 3)
ax2.set_ylim(-3, 3)
ax2.set_xlabel('acceleration (standard units)')
ax2.set_ylabel('msrp (standard units)')
ax2.grid(True, linestyle='--', alpha=0.7)
ax2.set_title('Acceleration vs MSRP')

plt.tight_layout()
plt.show()

结果如下：

这样我们就能更加清晰地观察到变量之间的关联性。由此可知，对于SUV车型，其售价与油耗呈负相关，而与加速度呈正相关。【不难理解，加速快的车性能更强，油耗更高，价格更贵】

回归线（Regression Line）

为简化讨论，下面我们只考虑标准化后的数据，并先考虑呈椭圆状分布的数据。示例：

import matplotlib.pyplot as plt
import numpy as np

np.random.seed(42) 
mean = [0, 0]
cov = [[1, 0.7], [0.7, 1]] 
x, y = np.random.multivariate_normal(mean, cov, 20000).T

fig, ax = plt.subplots(figsize=(6, 6),dpi=300)
ax.scatter(x, y, s=1, color='#00a8ff', alpha=0.8, marker='o')
ax.plot([-4, 4], [-4, 4], color='red', linewidth=2)
ax.plot([-4, 4], [-2.8, 2.8], color='green', linewidth=1.5)
ax.set_xlabel('x in standard units', fontsize=12)
ax.set_ylabel('y in standard units', fontsize=12)
ax.set_xlim(-4, 4)
ax.set_ylim(-4, 4)
ax.grid(True, linestyle='-', color='lightgrey')

plt.show()

correlation
可以看到，图中的绿线相比红线（ $y=x$ ）更能反映数据的均值分布。

实际上，绿线就被称为回归线，且其斜率与 $r$ 相等（即线的方程为 $y=rx$ ）。其反映的就是某一变量的平均水平与另一变量平均水平的线性关系。
于是我们也能反推出标准化前数据的回归线方程： $y - \overline{Y} = r \cdot \frac{s_Y}{s_X} (x - \overline{X}) \Longrightarrow y = \underbrace{r \cdot \frac{s_Y}{s_X} \cdot x}_{\text{斜率}} + \underbrace{(\overline{Y} - r \cdot \frac{s_Y}{s_X} \cdot \overline{X})}_{\text{截距}}$ 这也可以看作对数据的线性拟合。
让我们将身高预测例子的回归线进行补全：

# 前面代码省略（见Chapter 5）
from scipy import stats  

slope, intercept, r_value, p_value, std_err = stats.linregress(heights['Parent Average'], heights['Child'])
x_vals = np.array([heights['Parent Average'].min(), heights['Parent Average'].max()])
y_vals = slope * x_vals + intercept
plt.plot(x_vals, y_vals, color='red', linewidth=2, label=f'Regression line (r={slope:.2f})')
plt.xlabel('Parent Average Height', fontsize=12)
plt.ylabel('Child Height', fontsize=12)
plt.legend()
plt.grid(True, linestyle='--', alpha=0.6)
plt.show()

# install.packages("ggpmisc")
library(ggpmisc)  
heights$Prediction <- sapply(heights$Parent.Average, predict_child)
ggplot(heights, aes(x = Parent.Average)) +
  geom_point(aes(y = Child), alpha = 0.6, size = 3, shape = 21, 
             fill = "steelblue", colour = "white") +
  geom_point(aes(y = Prediction), alpha = 0.6, size = 3, shape = 21, 
             fill = "darkorange", colour = "white") +
  labs(x = "Parent Average Height", y = "Child Height") +
  geom_smooth(aes(y = Child), method = "lm", se = FALSE, 
              color = "red", size = 1.2) +
  stat_poly_eq(
    aes(y = Child,label = paste(..eq.label.., ..rr.label.., sep = "~~~~")),
    formula = y ~ x, parse = TRUE,
    label.x = "left", label.y = "top",  
    size = 4, color = "black"
  ) +
  theme_minimal()

结果如下：
regression-python regression-r
可以看到回归线与原来的均值预测非常接近，这说明线性回归的拟合效果很好。

最小二乘法（The Method of Least Squares）

现实中数据分布往往不是呈椭圆状的，那么此时如何用回归线进行拟合呢？此时我们就需要考虑对直线对样本的整体误差进行衡量。
从直觉上，我们会首先考虑预测值与实际值之差。但这样差值有正有负，难以聚合考量。于是我们就考虑另一种误差衡量方法——计算所有样本误差平方的均值，即均方误差（MSE）。
- 当然，为了与变量的单位保持一致，有时也对均方误差取平方根，即均方根误差（RMSE）。

可以证明，回归线对应的均方（根）误差是所有拟合直线中最小的。这里不作具体推导（应该会在《线性模型Cheat Sheet》中呈现），而直接使用程序进行数值优化：

import pandas as pd
import numpy as np
demographics = pd.read_csv('district_demographics2016.csv').drop(
    columns=['State', 'District', 'Percent voting for Clinton'])

def standard_units(x):
    """Converts an array x to standard units"""
    return (x - np.mean(x)) / np.std(x)

def correlation(t, x, y):
    x_su = standard_units(t[x])
    y_su = standard_units(t[y])
    return np.mean(x_su * y_su)

def slope(t, x, y):
    return correlation(t,x,y)*np.std(t[y])/np.std(t[x])
    
def intercept(t, x, y):
    return np.mean(t[y])-slope(t,x,y)*np.mean(t[x])

regression_slope = slope(demographics, 'College%', 'Median Income')
regression_intercept = intercept(demographics, 'College%', 'Median Income')
print(regression_slope, regression_intercept) # 回归线系数：1270.70168946 20802.5777667

def demographics_rmse(any_slope, any_intercept):
    x = demographics['College%']
    y = demographics['Median Income']
    estimate = any_slope*x + any_intercept
    return (np.mean((y - estimate) ** 2)) ** 0.5

print(minimize(demographics_rmse)) # 数值优化系数：1270.70168805  20802.57933807

可以看到，数值优化结果与回归线系数基本相同。因此，回归线有时也被称为“最小二乘线”。【这对于任意散点图分布都是成立的】

当然，利用minimize()函数，我们也可以进行非线性拟合，这里就不多赘述。另外，关于这个函数优化的原理，可以参考CS127或最优化相关知识。

视觉诊断

在得到回归线后，我们可以使用残差（即预测值与实际值之差）的分布评估预测的准确性。具体而言，我们可以通过绘制残差图评估线性回归的质量：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('default')
demographics = pd.read_csv('district_demographics2016.csv').drop(columns=['State', 'District', 'Percent voting for Clinton'])

def standard_units(x):
    """Converts an array x to standard units"""
    return (x - np.mean(x)) / np.std(x)

def correlation(t, x, y):
    x_su = standard_units(t[x])
    y_su = standard_units(t[y])
    return np.mean(x_su * y_su)

def slope(t, x, y):
    return correlation(t,x,y)*np.std(t[y])/np.std(t[x])
        
def intercept(t, x, y):
    return np.mean(t[y])-slope(t,x,y)*np.mean(t[x])

def fitted_values(t, x, y):
    """Return an array of the regression estimates (predictions) at all the x values"""
    a = slope(t, x, y)
    b = intercept(t, x, y)
    return a * t[x] + b

def residuals(t, x, y):
    predictions = fitted_values(t, x, y)
    return t[y] - predictions

demographics['Fitted Value'] = fitted_values(demographics, 'College%', 'Median Income')
demographics['Residual'] = residuals(demographics, 'College%', 'Median Income')

def plot_residuals(t, x, y):
    tbl = t.assign(
        Fitted=fitted_values(t, x, y),
        Residual=residuals(t, x, y)
    )

    # 图1：实际值 vs x 与 拟合值 vs x
    plt.figure(figsize=(10, 4),dpi=300)
    plt.subplot(1, 2, 1)
    plt.scatter(tbl[x], tbl[y], label='Actual', alpha=0.7)
    plt.scatter(tbl[x], tbl['Fitted'], label='Fitted', alpha=0.7)
    plt.xlabel(x)
    plt.ylabel(y)
    plt.legend()
    plt.title('Actual and Fitted Values')
    
    # 图2：残差 vs x
    plt.subplot(1, 2, 2)
    plt.scatter(tbl[x], tbl['Residual'], alpha=0.7)
    plt.axhline(y=0, color='r', linestyle='--')
    plt.xlabel(x)
    plt.ylabel('Residual')
    plt.title('Residuals vs. x')
    
    plt.tight_layout()
    plt.show()

plot_residuals(demographics, 'College%', 'Median Income')

residual-graph
可以看到，残差图中数据点在 $0$ 上下分布比较均匀，这说明回归模型的准确性在预测变量的不同取值范围内大致相同，因此这是一个好的回归。

当然，在更多情况下残差图用于诊断非线性关系，即当残差图中残差分布不均匀时，即认为数据之间可能有非线性关系。
除此之外，当数据在不同预测变量取值下分布方差不同（即准确性存在差异）时，称数据具有异方差性（Heteroscedasticity）。

数值诊断

除了通过残差图进行视觉诊断，我们还可以使用残差分布本身进行进一步考察：
- 首先，无论原始散点图形式如何，残差本身与预测变量都是没有线性相关性的，即二者相关系数为 $0$ ；
- 其次，残差的均值一定也为 $0$ ，而其标准差为 $\sqrt{1-r^2}\cdot s_Y^2$ 。由此可知，当相关系数 $r$ 越大，那么残差的标准差就越小，回归估计的误差就越低。
由此可以进一步得到以下结论：拟合预测值的标准差与 $s_Y^2$ 之比恰为 $|r|$ 。

解释

事实上，上述结论基于以下推导：约定符号

符号	含义
$SST$	因变量总离差平方和（即上面的 $s_Y^2$ ）
$SSR$	回归平方和（即拟合预测值的变异）
$SSE$	残差平方和（即残差的变异）
$Y_i$	第 $i$ 个样本点的因变量实际观测值
$\hat{Y}_i$	第 $i$ 个样本点的模型预测值（拟合值）
$\overline{Y}$	因变量实际值的样本均值

那么有

\begin{aligned} SST &= \sum (Y_i - \overline{Y})^2 \\ SSR &= \sum (\hat{Y}_i - \overline{Y})^2 \\ SSE &= \sum (Y_i - \hat{Y}_i)^2 \end{aligned}

且 $SST=SSR+SSE$ 。又因为 $\dfrac{SSE}{SST}=1-r^2$ ，因此 $\dfrac{SSR}{SSE}=r^2$ ，推得上述结论。

回归模型（Regression Model）

接下来，我们尝试将样本回归分析的结果推广到总体上。这涉及到对数据的假设，其集合即可称作模型。特别地，当假设数据点围绕某条直线随机波动，且这种波动可以用某种概率分布来描述时，那么这组假设就叫回归模型。
- 上述波动（即残差）分布要求为均值为 $0$ ，一般默认为正态分布（这样便于统计推断）。
如果总体的真实分布是一条直线上的值加上随机误差（噪声），那么当样本越大时，得到的回归线就与真实直线越接近，所以回归线可看作真实直线的良好近似。
- 那么当样本量受限时，应该如何评估回归线斜率与真实直线斜率的误差呢？是的，依旧可以使用Bootstrap自助法。
  - 具体而言，可以对散点图进行重采样，对每次重采样样本计算回归线斜率，最终得到Bootstrap分布。【后续就可以套用置信区间估计或假设检验】
  - 进一步，如果Bootstrap置信区间覆盖了 $0$ 点，那么数据之间就可能不存在线性关系，需要重新审视数据。
  - 在利用回归模型进行预测时，除了给出预测值之外，也可以使用Bootstrap得到预测值的置信区间。【注：预测变量值与均值越接近，其对应的Bootstrap置信区间就越窄】
总而言之，使用回归模型之前需要预先假设数据大致呈线性分布且分布均匀，这可以通过前面的视觉诊断和数值诊断实现。

分类（Classification）

略（主要讲的就是K近邻算法和一点点多元线性回归，总体较基础，笔者考虑在《机器学习方法》系列中进行补充）【见分类与聚类方法】

更新预测

同样略，这一块在《机器学习方法》中的朴素贝叶斯法中已经有所体现。（~~笔者想起课内上的《数据思维与实践》课程中也涉及了这一部分。~~）

ok，那么DATA 8系列就算完结了。可以看到，实际上笔者并没有像CS61A一样完整地过完所有内容，只是挑选了其中一些对笔者有用的部分，主要还是因为别的内容笔者或多或少都已经在其他课程中接触过。
总体上，个人感觉这门课更像是一门通识课，如果从这个角度上来看那么它确实比笔者学校里的那种通识课（或者计算机大课）好些。
另一方面，这门课也让我回想起去年暑假在UBC上的暑期班，果然国外课程对统计与数据科学的教授思路和国内很不一样呢（笑

原本笔者计划完结DATA 8后开坑DATA 100乃至DATA 140的，不过考虑到后续校内专业课涉及到的内容与其大概率重复，故决定DS方向暂时搁置，先更新CS相关课程。
以上~

DATA8 CHAPTER 10