线性代数/主题：最佳拟合直线

科学家经常会遇到一个没有解的系统，他们必须无论如何找到一个答案。也就是说，他们必须找到一个尽可能接近答案的值。

例如，假设我们有一枚硬币用来抛硬币。这枚硬币有某个比例$m$ 正面朝上的概率，由其物理构造决定，我们想知道$m$ 是否接近$1/2$。我们可以通过多次抛硬币来获得实验数据。这是对一枚美分硬币实验的结果，包括一些中间数字。

抛掷次数	30	60	90
正面朝上的次数	16	34	51

由于随机性，我们没有从这个样本中找到确切的比例 - 这个系统没有解。

${\begin{array}{*{1}{rc}r}30m&=&16\\60m&=&34\\90m&=&51\end{array}}$

也就是说，实验数据向量不在解的子空间中。

${\begin{pmatrix}16\\34\\51\end{pmatrix}}\not \in \{m{\begin{pmatrix}30\\60\\90\end{pmatrix}}\,{\big |}\,m\in \mathbb {R} \}$

但是，如上所述，我们希望找到$m$ 最接近有效的值。数据向量在直线子空间上的正交投影给出了我们的最佳猜测。

${\frac {{\begin{pmatrix}16\\34\\51\end{pmatrix}}\cdot {\begin{pmatrix}30\\60\\90\end{pmatrix}}}{{\begin{pmatrix}30\\60\\90\end{pmatrix}}\cdot {\begin{pmatrix}30\\60\\90\end{pmatrix}}}}\cdot {\begin{pmatrix}30\\60\\90\end{pmatrix}}={\frac {7110}{12600}}\cdot {\begin{pmatrix}30\\60\\90\end{pmatrix}}$

估计值（$m=7110/12600\approx 0.56$）有点高，但不是很高，所以这枚美分硬币可能足够公平。

斜率为$m\approx 0.56$ 的直线被称为此数据的最佳拟合直线。

最小化给定向量与用作右侧向量的向量之间的距离，最小化这些垂直长度的总和，因此我们说这条直线是通过最小二乘拟合得到的。

（这里垂直比例被放大了十倍，以使长度可见）。

我们整理了上面的方程，使得直线必须通过$(0,0)$，因为我们认为它是（我们对）斜率等于这枚硬币正面朝上概率的直线。我们也可以处理直线不必经过原点的案例。

例如，不同面额的美国货币有不同的平均流通时间（2美元钞票作为特殊情况被省略）。我们应该预期一张 25 美元的钞票能流通多久？

面额	1	5	10	20	50	100
平均寿命（年）	1.5	2	3	5	9	20

该图（见下图）看起来大致呈线性关系。它不是一条完美的直线，即线性方程组 $b+1m=1.5$，…，$b+100m=20$ 无解，但我们仍然可以使用正交投影来找到最佳近似解。考虑该线性方程组的系数矩阵及其常数向量，即实验确定的值。

$A={\begin{pmatrix}1&1\\1&5\\1&10\\1&20\\1&50\\1&100\end{pmatrix}}\qquad {\vec {v}}={\begin{pmatrix}1.5\\2\\3\\5\\9\\20\end{pmatrix}}$

投影到子空间部分的最后结果表明，系数 $b$ 和 $m$ 使得 $A$ 的列向量的线性组合尽可能接近向量 ${\vec {v}}$ 是 $({{A}^{\rm {trans}}}A)^{-1}{{A}^{\rm {trans}}}\cdot {\vec {v}}$ 的元素。一些计算得到截距为 $b=1.05$，斜率为 $m=0.18$。

将 $x=25$ 代入直线方程表明，这样的账单应该持续五年到六年之间。

最后，我们考虑男子一英里赛跑的时间 (Oakley & Baker 1977)。这些是给定年份1月1日生效的世界纪录。我们想预测何时会跑出3分40秒的成绩。

年份	1870	1880	1890	1900	1910	1920	1930
秒	268.8	264.5	258.4	255.6	255.6	252.6	250.4
年份	1940	1950	1960	1970	1980	1990	2000
秒	246.4	241.4	234.5	231.1	229.0	226.3	223.1

从下图可以看出，数据非常线性。有了这个输入

$A={\begin{pmatrix}1&1860\\1&1870\\\vdots &\vdots \\1&1990\\1&2000\end{pmatrix}}\qquad {\vec {v}}={\begin{pmatrix}280.0\\268.8\\\vdots \\226.3\\223.1\end{pmatrix}}$

本主题末尾的 Python 程序给出

${\text{slope}}=-0.35$ 和 ${\text{intercept}}=925.53$（四舍五入到两位小数；原始数据只精确到大约四分之一秒，因为其中大部分是手工计时）。

什么时候会跑完第 $220$ 英里的里程？求解最佳拟合线的方程可以估计出 $2008$ 年。

这个例子很有趣，但它也起到了警示作用——显然，数据的线性关系总有一天会失效（事实上，它在 1860 年之前就失效了）。

练习

这里的计算最好在计算机上完成。此外，一些问题需要更多的数据，这些数据可以在您的图书馆、网络、练习答案或练习后面的部分找到。

问题 1

使用最小二乘法判断这个实验中的硬币是否公平。

抛掷次数	8	16	24	32	40
正面次数	4	9	13	17	20

问题 2

对于男子一英里记录，我们并没有给出所有记录及其确切日期，而是通过周期性采样对数据进行了“平滑”。进行更长的计算，并将结论进行比较。

问题 3

找到男子 $1500$ 米跑的最佳拟合线。斜率与男子一英里的斜率相比如何？（距离很接近；一英里约为 $1609$ 米。）

问题 4

找到女子一英里记录的最佳拟合线。

问题 5

男子和女子一英里的最佳拟合线是否相交？

问题 6

1986 年挑战者号航天飞机爆炸时，人们对美国国家航空航天局决定发射的批评之一是，他们分析了密封圈失效次数与温度的关系（当然，密封圈失效导致了爆炸）。四个密封圈失效会导致火箭爆炸。美国国家航空航天局有 24 次之前飞行的数据。

温度（°F）	53	75	57	58	63	70	70	66	67	67	67
失效次数	3	2	1	1	1	1	1	0	0	0	0
温度（°F）	68	69	70	70	72	73	75	76	76	78	79	80	81
失效次数	0	0	0	0	0	0	0	0	0	0	0	0	0

当天的气温预报为 $31^{\circ }{\text{F}}$。

1. 美国国家航空航天局部分基于一张图表做出了发射决定，该图表只显示了至少发生一次密封圈失效的航班。找到最佳拟合这七次航班的线。根据这些数据，预测当温度为 $31$ 时，密封圈失效的次数，以及失效次数何时会超过四次。
2. 找到最佳拟合所有 24 次航班的线。根据这些额外的数据，预测当温度为 $31$ 时，密封圈失效的次数，以及失效次数何时会超过四次。

您认为哪种预测方法更准确？（在 (Dalal, Folkes & Hoadley 1989) 中有关于此的精彩讨论。）

问题 7

此表列出了前七颗行星到太阳的平均距离，以地球的平均距离为单位。

水星	金星	地球	火星	木星	土星	天王星
0.39	0.72	1.00	1.52	5.20	9.54	19.2

1. 绘制行星编号（水星为 $1$，等等）与距离的关系图。请注意，它不像一条直线，因此寻找最佳拟合线没有意义。
2. 然而，它看起来像一条指数曲线。因此，绘制行星编号与距离的对数的关系图。这看起来像一条直线吗？
3. 人们认为火星和木星之间的小行星带是曾经一颗行星破裂后留下的残骸。重新编号，使木星为 $6$，土星为 $7$，天王星为 $8$，并再次绘制与对数的关系图。这样看起来更好吗？
4. 对这些数据使用最小二乘法预测海王星的位置。
5. 重复上述步骤，预测冥王星的位置。
6. 这个公式对海王星和冥王星来说准确吗？

这种方法曾被用来帮助发现海王星（尽管第二条信息对历史的描述有误导性；实际上，海王星在位置 $9$ 被发现，促使人们在位置 $5$ 寻找“失踪的行星”。参见 (Gardner 1970)

问题 8

威廉·贝内特提出了美国领先文化指标指数 (Bennett 1993)。引用的统计数据包括平均每天看电视的时间和平均综合 SAT 成绩。

	1960	1965	1970	1975	1980	1985	1990	1992
电视	5:06	5:29	5:56	6:07	6:36	7:07	6:55	7:04
SAT	975	969	948	910	890	906	900	899

假设在观看电视的时间和 SAT 成绩下降之间存在因果关系（在这篇文章中，贝内特先生并没有主张两者之间存在直接联系）。

1. 找到将平均每日观看电视时间（自变量）与 SAT 成绩（因变量）联系起来的最佳拟合线。
2. 找到最近的平均每日观看电视时间估计值（贝内特引用尼尔森媒体研究作为这些估计值的来源）。估计相关的 SAT 成绩。你的估计值与实际平均值有多接近？（警告：SAT 最近进行了更改，因此你应该调查是否需要对报告的平均值进行一些调整才能进行有效的比较）。

解决方案

计算机代码

#!/usr/bin/python
# least_squares.py   calculate the line of best fit for a data set
# data file format: each line is two numbers, x and y
n = 0
sum_x = 0
sum_y = 0
sum_x_squared = 0
sum_xy = 0

fn = raw_input("Name of the data file? ")
datafile = open(fn,"r")
while 1:
  ln = datafile.readline()
  if ln:
    data = ln.split()
    x = float(data[0])
    y = float(data[1])
    n += 1
    sum_x += x
    sum_y += y
    sum_x_squared += x*x
    sum_xy += x*y
  else:
    break
datafile.close()

slope = (n*sum_xy - sum_x*sum_y) / (n*sum_x_squared - sum_x**2)
intercept = (sum_y - slope*sum_x)/n
print "line of best fit: slope= %f  intercept= %f" % (slope, intercept)

附加数据

以下是一些关于世界纪录发展的数据（取自跑者世界网站）。

男子英里纪录的发展
时间	姓名	日期
4:52.0	Cadet Marshall (GBR)	02Sep52
4:45.0	Thomas Finch (GBR)	03Nov58
4:40.0	Gerald Surman (GBR)	24Nov59
4:33.0	George Farran (IRL)	23May62
4:29 3/5	Walter Chinnery (GBR)	10Mar68
4:28 4/5	William Gibbs (GBR)	03Apr68
4:28 3/5	Charles Gunton (GBR)	31Mar73
4:26.0	Walter Slade (GBR)	30May74
4:24 1/2	Walter Slade (GBR)	19Jun75
4:23 1/5	Walter George (GBR)	16Aug80
4:19 2/5	Walter George (GBR)	03Jun82
4:18 2/5	Walter George (GBR)	21Jun84
4:17 4/5	Thomas Conneff (USA)	26Aug93
4:17.0	Fred Bacon (GBR)	06Jul95
4:15 3/5	Thomas Conneff (USA)	28Aug95
4:15 2/5	John Paul Jones (USA)	27May11
4:14.4	John Paul Jones (USA)	31May13
4:12.6	Norman Taber (USA)	16Jul15
4:10.4	Paavo Nurmi (FIN)	23Aug23
4:09 1/5	Jules Ladoumegue (FRA)	04Oct31
4:07.6	Jack Lovelock (NZL)	15Jul33
4:06.8	Glenn Cunningham (USA)	16Jun34
4:06.4	Sydney Wooderson (GBR)	28Aug37
4:06.2	Gunder Hagg (SWE)	01Jul42
4:04.6	Gunder Hagg (SWE)	04Sep42
4:02.6	Arne Andersson (SWE)	01Jul43
4:01.6	Arne Andersson (SWE)	18Jul44
4:01.4	Gunder Hagg (SWE)	17Jul45
3:59.4	Roger Bannister (GBR)	06May54
3:58.0	John Landy (AUS)	21Jun54
3:57.2	Derek Ibbotson (GBR)	19Jul57
3:54.5	Herb Elliott (AUS)	06Aug58
3:54.4	Peter Snell (NZL)	27Jan62
3:54.1	Peter Snell (NZL)	17Nov64
3:53.6	Michel Jazy (FRA)	09Jun65
3:51.3	Jim Ryun (USA)	17Jul66
3:51.1	Jim Ryun (USA)	23Jun67
3:51.0	Filbert Bayi (TAN)	17May75
3:49.4	John Walker (NZL)	12Aug75
3:49.0	Sebastian Coe (GBR)	17Jul79
3:48.8	Steve Ovett (GBR)	01Jul80
3:48.53	Sebastian Coe (GBR)	19Aug81
3:48.40	Steve Ovett (GBR)	26Aug81
3:47.33	Sebastian Coe (GBR)	28Aug81
3:46.32	Steve Cram (GBR)	27Jul85
3:44.39	Noureddine Morceli (ALG)	05Sep93
3:43.13	Hicham el Guerrouj (MOR)	07Jul99

男子 1500 米纪录的发展
时间	姓名	日期
4:09.0	John Bray (USA)	30May00
4:06.2	Charles Bennett (GBR)	15Jul00
4:05.4	James Lightbody (USA)	03Sep04
3:59.8	Harold Wilson (GBR)	30May08
3:59.2	Abel Kiviat (USA)	26May12
3:56.8	Abel Kiviat (USA)	02Jun12
3:55.8	Abel Kiviat (USA)	08Jun12
3:55.0	Norman Taber (USA)	16Jul15
3:54.7	John Zander (SWE)	05Aug17
3:53.0	Paavo Nurmi (FIN)	23Aug23
3:52.6	Paavo Nurmi (FIN)	19Jun24
3:51.0	Otto Peltzer (GER)	11Sep26
3:49.2	Jules Ladoumegue (FRA)	05Oct30
3:49.0	Luigi Beccali (ITA)	17Sep33
3:48.8	William Bonthron (USA)	30Jun34
3:47.8	Jack Lovelock (NZL)	06Aug36
3:47.6	Gunder Hagg (SWE)	10Aug41
3:45.8	Gunder Hagg (SWE)	17Jul42
3:45.0	Arne Andersson (SWE)	17Aug43
3:43.0	Gunder Hagg (SWE)	07Jul44
3:42.8	Wes Santee (USA)	04Jun54
3:41.8	John Landy (AUS)	21Jun54
3:40.8	Sandor Iharos (HUN)	28Jul55
3:40.6	Istvan Rozsavolgyi (HUN)	03Aug56
3:40.2	Olavi Salsola (FIN)	11Jul57
3:38.1	Stanislav Jungwirth (CZE)	12Jul57
3:36.0	Herb Elliott (AUS)	28Aug58
3:35.6	Herb Elliott (AUS)	06Sep60
3:33.1	Jim Ryun (USA)	08Jul67
3:32.2	Filbert Bayi (TAN)	02Feb74
3:32.1	Sebastian Coe (GBR)	15Aug79
3:31.36	Steve Ovett (GBR)	27Aug80
3:31.24	Sydney Maree (usa)	28Aug83
3:30.77	Steve Ovett (GBR)	04Sep83
3:29.67	Steve Cram (GBR)	16Jul85
3:29.46	Said Aouita (MOR)	23Aug85
3:28.86	Noureddine Morceli (ALG)	06Sep92
3:27.37	Noureddine Morceli (ALG)	12Jul95
3:26.00	Hicham el Guerrouj (MOR)	14Jul98

女子英里纪录的发展
时间	姓名	日期
6:13.2	Elizabeth Atkinson (GBR)	24Jun21
5:27.5	Ruth Christmas (GBR)	20Aug32
5:24.0	Gladys Lunn (GBR)	01Jun36
5:23.0	Gladys Lunn (GBR)	18Jul36
5:20.8	Gladys Lunn (GBR)	08May37
5:17.0	Gladys Lunn (GBR)	07Aug37
5:15.3	Evelyne Forster (GBR)	22Jul39
5:11.0	Anne Oliver (GBR)	14Jun52
5:09.8	Enid Harding (GBR)	04Jul53
5:08.0	Anne Oliver (GBR)	12Sep53
5:02.6	Diane Leather (GBR)	30Sep53
5:00.3	Edith Treybal (ROM)	01Nov53
5:00.2	Diane Leather (GBR)	26May54
4:59.6	Diane Leather (GBR)	29May54
4:50.8	Diane Leather (GBR)	24May55
4:45.0	Diane Leather (GBR)	21Sep55
4:41.4	Marise Chamberlain (NZL)	08Dec62
4:39.2	Anne Smith (GBR)	13May67
4:37.0	Anne Smith (GBR)	03Jun67
4:36.8	Maria Gommers (HOL)	14Jun69
4:35.3	Ellen Tittel (FRG)	20Aug71
4:34.9	Glenda Reiser (CAN)	07Jul73
4:29.5	Paola Pigni-Cacchi (ITA)	08Aug73
4:23.8	Natalia Marasescu (ROM)	21May77
4:22.1	Natalia Marasescu (ROM)	27Jan79
4:21.7	玛丽·德克尔（美国）	1980年1月26日
4:20.89	卢德米拉·维塞尔科娃（苏联）	1981年9月12日
4:18.08	玛丽·德克尔-塔布（美国）	1982年7月9日
4:17.44	马里卡·普伊卡（罗马尼亚）	1982年9月16日
4:15.8	纳塔莉亚·阿特约莫娃（苏联）	1984年8月5日
4:16.71	玛丽·德克尔-斯莱尼（美国）	1985年8月21日
4:15.61	保拉·伊万（罗马尼亚）	1989年7月10日
4:12.56	斯维特拉娜·马斯特科娃（俄罗斯）	1996年8月14日

参考文献

• Bennett, William (1993年3月15日), "量化美国衰退", 华尔街日报{{citation}}: CS1 maint: date and year (link)
• Dalal, Siddhartha; Folkes, Edward; Hoadley, Bruce (1989年秋季), "从挑战者号中学到的教训：统计学视角", 统计学学生杂志, pp. 14–18{{citation}}: CS1 maint: date and year (link)
• Gardner, Martin (1970年4月), "数学游戏，太阳系中嵌入的一些数学奇点", 科学美国人, pp. 108–112{{citation}}: CS1 maint: date and year (link)
• Oakley, Cletus; Baker, Justine (1977年4月), "最小二乘法与3:40英里", 数学教师{{citation}}: CS1 maint: date and year (link)