最接近直線

作者：王一哲
日期：2019/8/5

前言

物理科及化學科在處理實驗數據時，通常會先繪製 XY 散佈圖，如果資料點看起來有線性的關係，就要用最小平方法計算最接近直線的斜率和截距，再試著解釋斜率和截距的物理意義。但如果更講究一點，由於測量數據時本身就會有一定的不準量，最接近直線的斜率、截距也會有不準量，要用什麼方法才能很方便地算出這些數呢？

最接近直線公式

以下的公式是我從物理奧林匹亞培訓講義上抄來的。假設操縱變因（自變數）為$x$、不準量為$\Delta x$，測量的應變變因（應變數）為$y$、不準量為$\Delta y$，數據共有$n$組。則最接近直線

$$斜率 \quad a = \frac{\sum x \sum y - n \sum xy}{(\sum x)^2 - n \sum x^2}$$
$$截距 \quad b = \frac{\sum x \sum xy - \sum y \sum x^2}{(\sum x)^2 - n \sum x^2}$$
$$相關係數 \quad R = \frac{\sum x \sum y - n \sum xy}{\sqrt{\left[ (\sum x)^2 - n \sum x^2 \right] \left[ (\sum y)^2 - n \sum y^2 \right]}}$$

為了要計算斜率及截距的不準量，我們先定義下列的符號

$$\sigma_x = \sqrt{\frac{\sum (\Delta x)^2}{n}}$$
$$\sigma_y = \sqrt{\frac{\sum (\Delta y)^2}{n}}$$
$$\sigma = \sqrt{\sigma_y^2 + a^2 \sigma_x^2}$$

斜率與截距的不準量分別為

$$斜率的不準量 \quad \Delta a = \sqrt{\frac{n \sigma^2}{n \sum x^2 - (\sum x)^2}}$$
$$截距的不準量 \quad \Delta b = \sqrt{\frac{\sigma^2 \sum x^2}{n \sum x^2 - (\sum x)^2}}$$

如果要將$x$軸的不準量$\Delta x$加到$y$軸上，等效的$y$軸不準量計算方式為

$$\Delta y_{eff} = \sqrt{(\Delta x)^2 + (a \Delta y)^2}$$

使用 LibreOffice Calc

計算最接近直線斜率、截距及不準量

首先我們用 Calc 產生一些資料，假設 $y=ax+b$、$a=2$、$b=3$，再加上亂數使資料有一點偏差，產生的資料表如下：

index	x	$\Delta x$	y	$\Delta y$
0	1.0	0.1	4.9	0.5
1	2.0	0.1	6.9	0.5
2	3.0	0.1	9.0	0.5
3	4.0	0.1	11.1	0.6
4	5.0	0.1	13.2	0.6
5	6.0	0.2	15.2	0.8
6	7.0	0.2	16.9	0.8
7	8.0	0.2	18.9	0.8
8	9.0	0.2	20.9	0.8
9	10.0	0.2	22.8	0.8

如果不想要手動輸入資料，可以直接將網頁上的資料選取並複製，直接貼到試算表中。接著我們先計算一些需要用到的數值，例如$x^2$、$y^2$、$xy$、$(\Delta x)^2$、$(\Delta y)^2$，作法相當類似，以$x^2$為例，先用滑鼠左鍵點選儲存格F2，接著輸入

=B2^2

再將滑鼠游標移到儲存格F2的右下角，當滑鼠游標變為黑色十字時按住滑鼠左鍵向下拖曳到儲存格F11，Calc會自動以公式B3^2、B4^2……這樣的公式計算每個儲存格對應的數值。

接著於儲存格A14到儲存格J14計算對應的數量及總合。$n$是資料點數量，先用滑鼠左鍵點選儲存格A14，接著輸入

=COUNT(A2:A11)

$\sum x$是$x$數值總合，先用滑鼠左鍵點選儲存格B14，接著輸入

=SUM(B2:B11)

以上兩個函式的詳細用法請參考官方說明書：COUNT、SUM。

計算需要使用的數值

我們可以用 Calc 內建的函式計算最接近直線的斜率，先用滑鼠左鍵點選儲存格M2，接著輸入

=SLOPE(D2:D11,B2:B11)

語法為

=SLOPE(y軸資料, x軸資料)

如果要手動輸入公式，先用滑鼠左鍵點選儲存格N2，接著輸入

=(B14*D14-A14*H14)/(B14^2-A14*F14)

兩者計算結果均為1.98787878787879。



再來計算最接近直線的截距，先用滑鼠左鍵點選儲存格M3，接著輸入

=INTERCEPT(D2:D11,B2:B11)

語法為

=INTERCEPT(y軸資料, x軸資料)

如果要手動輸入公式，先用滑鼠左鍵點選儲存格N3，接著輸入

=(B14*H14-D14*F14)/(B14^2-A14*F14)

兩者計算結果約為3.04666666666668，直到小數點後第14位才有差異。



最後計算最接近直線的$R^2$值，先用滑鼠左鍵點選儲存格M9，接著輸入

=RSQ(D2:D11,B2:B11)

語法為

=RSQ(y軸資料, x軸資料)

如果要手動輸入公式，先用滑鼠左鍵點選儲存格N9，接著輸入

=((B14*D14-A14*H14)/SQRT((B14^2-A14*F14)*(D14^2-A14*G14)))^2

兩者計算結果約為0.999497575579201，直到小數點後第15位才有差異。以上三個函式的詳細用法請參考官方說明書。

但是 Cacl 沒有內建計算斜率、截距不準量的函數，我們只能手動輸入公式。先計算$\sigma_x$、$\sigma_y$、$\sigma$，公式分別為

=SQRT(I14/A14)
=SQRT(J14/A14)
=SQRT(N5^2+N2^2*N4^2)

斜率不準量$\Delta a$公式為

=SQRT(A14*N6^2/(A14*F14-B14^2))

計算結果為0.082813521943926。截距不準量$\Delta b$公式為

=SQRT(N6^2*F14/(A14*F14-B14^2))

計算結果為0.513844390399612。

繪製帶有誤差槓的XY散佈圖

雖然現在正式的名稱已經改為不準量 (uncertainty) ，不過一般而言還是會稱為誤差 (error)，在繪製XY散佈圖時資料點要加上誤差槓 (error bar)。

插入XY散佈圖

首先選取儲存格D1:D11（下圖1）、B1:B11（下圖2），再點右上角的插入圖表（下圖3），或是選取資料後點選單中的插入⇒圖表，插入XY散佈圖。

由工具列插入圖表

由選單插入圖表

接著會跳出圖表精靈，選擇圖表類型中的XY(散佈)，種類為僅限點，再按下一步。

圖表精靈圖表類型

由於我們已經先選取了資料範圍，不需要再指定資料範圍。選擇以欄表示的為資料序列，勾選以第一列作為標籤，再按下一步。

圖表精靈資料範圍

由於我們目前只想在這張圖中畫出一組資料序列，直接按下一步即可。

圖表精靈資料序列

我們習慣上不會加上題名和副題，會標示在圖的說明文字中，這兩格空白即可。X軸、Y軸這兩格寫上代表的物理量及單位。由於圖上只有一組數據，所以不需要顯示圖例。格線的部分則是依照個人喜好決定，沒有一定的作法。按完成之後產生的圖形如下。

圖表精靈圖表元素

只有資料點的XY散佈圖

修改圖片格式

接下來還可以再修改資料點的格式，先在圖上連點按滑鼠左鍵兩下進入圖表，再於資料點上按滑鼠右鍵開啟快速選單，選取設定資料序列格式。

快速選單

在資料序列視窗中的線條分頁，可以修改線條的樣式及顏色，也可以修改資料點圖示。

資料序列視窗

從右側的選取⇒圖庫，選取自己喜歡的樣式。

資料序列圖示

修改資料點圖示的XY散佈圖

如果覺得預設的座標軸標示、標籤的字體太小，可以在軸上面按滑鼠右鍵，從選單中選取設定軸格式。

設定軸格式

選擇自己喜歡的字型、樣式及大小，再按下確定即可。

設定Y軸字型

我們可以用類似的方法修改座標軸標籤。

修改座標軸標籤

設定Y軸標題字型

修改座標軸格式和標籤後的XY散佈圖

插入誤差槓

在資料點上按滑鼠右鍵，從選單中選取插入X誤差線。

插入X誤差線選單

由於我們的C欄資料就是$x$軸不準量，所以我們要選擇誤差類別中的儲存格範圍（下圖1），再按下圖2的標籤從資料表中選取儲存格C2:C11，由於我們的不準量正負值相同，只要勾選下圖3的兩者同值，最後按下確定（下圖4）即可。接著可以用類似的方法插入Y誤差線。

資料序列「y」的X誤差線

插入誤差線後的XY散佈圖

插入趨勢線

如果要在圖上加入趨勢線，可以在資料點上按滑鼠右鍵，從選單中選取插入趨勢線。在下圖的視窗中選取線性（下圖1），如果要在圖上顯示方程式及決定係數$R^2$，則要勾選圖中標籤2、3，最後再按下確定（下圖4）即可。但是我們通常不會在圖上顯示方程式及決定係數，會另外標示於圖表說明中。

資料序列「y」的趨勢線

插入擬合線後的XY散佈圖

匯出成圖檔

如果要把XY散佈圖匯出成圖檔，不要進入編輯圖表的模式，直接於圖上按滑鼠右鍵，從選單中選取匯出為影像，決定存檔路徑、檔名及圖片格式，最後按存檔。

匯出為影像選單

決定存檔路徑、檔名及圖片格式

匯出後的影像

使用 SciDAVis

從 CSV 檔匯入資料

先將以下的資料複製、貼上到文字編輯器中，將檔案儲存為data.csv，文字編碼應該用 ASCII 或 UTF-8 皆可。

index,x,errx,y,erry
0,1.0,0.1,4.9,0.5
1,2.0,0.1,6.9,0.5
2,3.0,0.1,9.0,0.5
3,4.0,0.1,11.1,0.6
4,5.0,0.1,13.2,0.6
5,6.0,0.2,15.2,0.8
6,7.0,0.2,16.9,0.8
7,8.0,0.2,18.9,0.8
8,9.0,0.2,20.9,0.8
9,10.0,0.2,22.8,0.8

由選單：File⇒Import ASCII 或 Ctrl+K，開啟匯入檔案視窗。

匯入檔案選單

由於資料是以逗號分隔，在匯入檔案視窗中，要記得將 Separator 選項改成逗號。另外也可以選擇要將資料匯入哪一個表格，可以蓋掉現有表格的資料，也可以另外增表格。

匯入檔案視窗

匯入資料檔後的表格

設定欄位資料格式

接著在要修改的欄位上按滑鼠左鍵，在畫面右側也可以更改欄位的描述及資料格式。

更改欄位的描述

更改欄位的資料格式

在要修改的欄位上按滑鼠右鍵，選擇Set Column(s) As即可修改欄位的資料用途，分別將欄位index設為X、欄位x設為X、欄位errx設為X Error、欄位y設為Y、欄位erry設為Y Error。

設定欄位的資料用途

繪製帶有誤差槓的XY散佈圖

插入XY散佈圖

先選取欄位x、y之後，再從選單選取Plot⇒Scatter，或是在欄位上按滑鼠右鍵，從快速選單中選取Plot⇒Scatter。

從選單插入XY散佈圖

從快速選單插入XY散佈圖

預設格式的XY散佈圖

修改圖片格式

如果想要修改資料點格式，可以在圖片的資料點上連點滑鼠左鍵兩下，勾選Fill Color會在點中填滿指定的顏色，Edge Color是點的邊緣線條顏色。

修改資料點格式

點選Plot Associations，可以修改資料點x軸、y軸數值對應的欄位。

修改資料點對應的欄

如果想要修改座標軸及標籤的格式，可以在圖片的座標軸上連點滑鼠左鍵兩下，點選Axis分頁，在Title欄位中可以修改標籤內容，支援HTML語法；點選Axis Font可以修改座標軸數字的字型及大小。

修改座標軸性質

修改座標軸數字的字型

修改格式後的XY散佈圖

插入誤差槓

點選圖片後從選單選取Graph⇒Add Error Bars或按快速鍵Ctrl+B。

插入誤差槓

勾選Existing column，從資料表中選取不準量的來源，在視窗的下方勾選X Error Bars或Y Error Bars，最後點選右上方的Add。

插入誤差槓的XY散佈圖

插入趨勢線

點選圖片後從選單選取Analysis⇒Quick Fit⇒Fit Linear。

選單插入趨勢線

擬合的結果如下，看起來只有考慮y軸不準量。

$$ a = 1.9979669965615 \pm 0.0715310360117212 $$
$$ b = 2.99177166800124 \pm 0.371842523931618 $$
$$ R^2 = 0.999509982857237 $$

擬合結果

如果要再開啟輸出結果的視窗，可以從選單選取View⇒Results Log。

開啟輸出結果視窗

匯出成圖檔

點選圖片後從選單選取File⇒Export Graph⇒Current或按快速鍵Alt+G。

匯出圖檔選單

決定存檔路徑、檔名及圖片格式，最後按存檔。經過測試後發現，存成png檔大約400KB，存成jpg檔大約35KB，但是在電腦螢幕上看起來效果差不多。

匯出圖檔視窗

匯出的XY散佈圖

使用 Python

最後我們用Python處理數據，需要用到的函式庫為NumPy、SciPy及Matplotlib。程式碼如下

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# 從 data.csv 讀取資料
data = np.loadtxt("data.csv", dtype=float, delimiter=',', skiprows=1)
x = data[:, 1]
errx = data[:, 2]
y = data[:, 3]
erry = data[:, 4]

"""
方法1: 手動輸入公式
"""
n = len(x)
sum_x = x.sum()
sum_y = y.sum()
sum_xy = np.sum(x*y)
sum_x2 = np.sum(x**2)
sum_y2 = np.sum(y**2)
sum_errx2 = np.sum(errx**2)
sum_erry2 = np.sum(erry**2)
sigma_x = np.sqrt(sum_errx2/n)
sigma_y = np.sqrt(sum_erry2/n)

a = (sum_x*sum_y - n*sum_xy) / (sum_x**2 - n*sum_x2)
b = (sum_x*sum_xy - sum_y*sum_x2) / (sum_x**2 - n*sum_x2)
sigma = np.sqrt(sigma_y**2 + a**2 * sigma_x**2)
erra = np.sqrt(n*sigma**2 / (n*sum_x2 - sum_x**2))
errb = np.sqrt(sigma**2 * sum_x2 / (n*sum_x2 - sum_x**2))
r = (sum_x*sum_y - n*sum_xy) / np.sqrt((sum_x**2 - n*sum_x2) * (sum_y**2 - n*sum_y2))

print("方法1: 手動輸入公式")
print("slope: %.12f +/- %.12f" % (a, erra))
print("intercept: %.12f +/- %.12f" % (b, errb))
print("r-squared: %.12f" % r**2)
print("sigma_x: %.12f   sigma_y: %.12f   sigma: %.12f" % (sigma_x, sigma_y, sigma))

"""
方法2: 使用內建函數
"""

# 計算最接近直線的斜率、截距及不準量
popt, pcov = curve_fit(lambda x, a, b: a*x+b, x, y, sigma=erry, absolute_sigma=True)

print("方法2: 使用內建函數")
print("slope: %.12f +/- %.12f" % (popt[0], np.sqrt(pcov[0, 0])))
print("intercept: %.12f +/- %.12f" % (popt[1], np.sqrt(pcov[1, 1])))

# 計算最接近直線的R平方值
cor = np.corrcoef(x, y)[0, 1]
print("r-squared: %.12f" % cor**2)

# 建立繪圖物件
plt.figure(figsize=(6, 4.5), dpi=72)
plt.xlabel('x', fontsize=14)
plt.ylabel('y', fontsize=14)
plt.grid(color='grey', linestyle='--', linewidth=1)
plt.plot(x, y, color='blue', marker='o', markersize=10, linestyle='')

# 加上 error bar
plt.errorbar(x, y, xerr=errx, yerr=erry, color='black', fmt='o', capsize=4, markersize=0, linewidth=1)

# 加上最接近直線
line = popt[0]*x + popt[1]
plt.plot(x, line, color='red', marker='', linestyle='-', linewidth=2)

# 儲存、顯示圖片
plt.savefig('linear.svg')
plt.savefig('linear.png')
plt.show()

從 CSV 檔匯入資料

使用的函式為numpy.loadtxt，詳細的用法請參考官方說明書。程式碼第6行

data = np.loadtxt("data.csv", dtype=float, delimiter=',', skiprows=1)

意思是從檔案data.csv讀取資料，資料以逗號分隔，加上skiprows=1略過第1列，將資料指定給物件data，資料格式為 NumPy Ndarray。程式碼第7到10行，再將data的第1、2、3、4欄的資料分別存入1維陣列x、errx、y、erry。

方法1: 手動輸入公式

我們先手動算出需要用到的數值。例如程式碼第15行，用函式len計算資料數量再指定物件n。語法有兩種

[陣列物件].len
len([陣列物件])

程式碼第16到22行，用函式numpy.sum計算陣列資料的總合，由於資料為1維陣列語法為

[陣列物件].sum
np.sum([陣列物件])

詳細的用法請參考官方說明書。



程式碼第26到30行，將數值代入文章最前面的公式，計算結果為

$$ a = 1.987878787879 \pm 0.082813521944 $$
$$ b = 3.046666666667 \pm 0.513844390400 $$
$$ R^2 = 0.999497575579 $$

方法2: 使用內建函數

我們使用scipy.optimize.curve_fit計算最接近直線的斜率、截距及不準量，由於程式一開始已經引入這個函式，所以語法為

curve_fit([擬合函數型式], [x軸資料], [y軸資料], sigma=[y軸不準量], absolute_sigma=True)

我們在程式中寫成

popt, pcov = curve_fit(lambda x, a, b: a*x+b, x, y, sigma=erry, absolute_sigma=True)

將輸出的資料拆成兩部分，分別指定給物件popt和pcov。利用lambda運算式定義函式a*x+b，可以少用一次def。選項absolute_sigma預設值為False，若設定為True在計算時會考慮y軸資料的不準量。詳細的用法請參考官方說明書。



接著使用numpy.corrcoef計算相關係數$R$，語法為

np.corrcoef([x軸資料], [y軸資料])

回傳值是一個2*2的陣列，我們要的$R^2$陣列右上角元素值的平方，因此程式碼寫為

cor = np.corrcoef(x, y)[0, 1]
print("r-squared: %.12f" % cor**2)

詳細的用法請參考官方說明書。



計算結果為

$$ a = 1.997966995944 \pm 0.071531035658 $$
$$ b = 2.991771672485 \pm 0.371842522314 $$
$$ R^2 = 0.999497575579 $$

繪製帶有誤差槓的XY散佈圖

最後我們用matplotlib.pyplot繪圖，由於函式庫的全名太長，通常引入時會寫成

import matplotlib.pyplot as plt

固定的簡稱為plt。程式碼的說明如下

1. 第55行：figure，開啟繪圖物件，寬度為6英吋、高度為4.5英吋、dpi為72。
2. 第56行：xlabel，設定x軸標籤內容為x、字體大小為14號。
3. 第57行：ylabel，設定y軸標籤內容為y、字體大小為14號。
4. 第58行：grid，設定格線為寬度1的灰色虛線。
5. 第59行：plot，以1維陣列x、y的資料繪圖，只有資料點，點的格式為藍色圓形，大小為10。
6. 第62行：errorbar，加上誤差槓，其中xerr為x軸不準量、yerr為x軸不準量，color為線條顏色，capsize為槓的寬度，linewidth為線條寬度，由於我不想要再畫一次點，因此將markersize設為0，詳細的用法請參考官方說明書。
7. 第65、66行：為了加上最接近直線，我們先用擬合得到的斜率popt[0]、截距popt[1]，將陣列x資料代入計算對應的y值，再用plot函式畫線。

用 matplotlib 繪製的XY散佈圖

結語

以下是不同工具的計算結果，其中SciDAVis和Python(內建)在擬合時會考慮y軸不準量，而且計算結果相當接近，如果考慮學生學習工具上的難度，我會建議學生使用SciDAVis，如果不排斥看到程式碼再學著使用Python。我相信按照這個標準得到的XY散佈圖，應該可以符合科展評審的標準。

	內建公式	手動輸入	SciDAVis	Python(手動)	Python(內建)
斜率a	1.98787878787879	1.98787878787879	1.99796699656150	1.987878787879	1.997966995944
截距b	3.04666666666667	3.04666666666668	2.99177166800124	3.046666666667	2.991771672485
$\sigma_x$		0.158113883008419		0.158113883008
$\sigma_y$		0.68337398253079		0.683373982531
$\sigma$		0.7521911671127649		0.752191167113
不準量$\Delta a$		0.082813521943926	0.071531036011721	0.082813521944	0.071531035658
不準量$\Delta b$		0.513844390399612	0.371842523931618	0.513844390400	0.371842522314
$R^2$	0.999497575579200	0.999497575579201	0.999509982857237	0.999497575579	0.999497575579

---

HackMD 版本連結：https://hackmd.io/@yizhewang/Sy6OHP2GS