日期：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
接著於儲存格A14到儲存格J14計算對應的數量及總合。$n$是資料點數量，先用滑鼠左鍵點選儲存格A14，接著輸入
=COUNT(A2:A11)
=SUM(B2:B11)
計算需要使用的數值
我們可以用 Calc 內建的函式計算最接近直線的斜率，先用滑鼠左鍵點選儲存格M2，接著輸入
=SLOPE(D2:D11,B2:B11)
=SLOPE(y軸資料, x軸資料)
=(B14*D14-A14*H14)/(B14^2-A14*F14)
再來計算最接近直線的截距，先用滑鼠左鍵點選儲存格M3，接著輸入
=INTERCEPT(D2:D11,B2:B11)
=INTERCEPT(y軸資料, x軸資料)
=(B14*H14-D14*F14)/(B14^2-A14*F14)
最後計算最接近直線的$R^2$值，先用滑鼠左鍵點選儲存格M9，接著輸入
=RSQ(D2:D11,B2:B11)
=RSQ(y軸資料, x軸資料)
=((B14*D14-A14*H14)/SQRT((B14^2-A14*F14)*(D14^2-A14*G14)))^2
但是 Cacl 沒有內建計算斜率、截距不準量的函數，我們只能手動輸入公式。先計算$\sigma_x$、$\sigma_y$、$\sigma$，公式分別為
=SQRT(I14/A14)
=SQRT(J14/A14)
=SQRT(N5^2+N2^2*N4^2)
=SQRT(A14*N6^2/(A14*F14-B14^2))
=SQRT(N6^2*F14/(A14*F14-B14^2))
繪製帶有誤差槓的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。程式碼的說明如下
- 第55行：figure，開啟繪圖物件，寬度為6英吋、高度為4.5英吋、dpi為72。
- 第56行：xlabel，設定x軸標籤內容為x、字體大小為14號。
- 第57行：ylabel，設定y軸標籤內容為y、字體大小為14號。
- 第58行：grid，設定格線為寬度1的灰色虛線。
- 第59行：plot，以1維陣列x、y的資料繪圖，只有資料點，點的格式為藍色圓形，大小為10。
- 第62行：errorbar，加上誤差槓，其中xerr為x軸不準量、yerr為x軸不準量，color為線條顏色，capsize為槓的寬度，linewidth為線條寬度，由於我不想要再畫一次點，因此將markersize設為0，詳細的用法請參考官方說明書。
- 第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
沒有留言:
張貼留言