線性代數/主題：最佳擬合直線

科學家經常會遇到一個沒有解的系統，他們必須無論如何找到一個答案。也就是說，他們必須找到一個儘可能接近答案的值。

例如，假設我們有一枚硬幣用來拋硬幣。這枚硬幣有某個比例$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)