<style>
.markdown-body table{
display: unset;
}
</style>
# 最接近直線
> 作者：王一哲
> 日期：2019/8/5
<br />
## 前言
物理科及化學科在處理實驗數據時，通常會先繪製 XY 散佈圖，如果資料點看起來有線性的關係，就要用最小平方法計算最接近直線的斜率和截距，再試著解釋斜率和截距的物理意義。但如果更講究一點，由於測量數據時本身就會有一定的不準量，最接近直線的斜率、截距也會有不準量，要用什麼方法才能很方便地算出這些數呢？
<br />
## 最接近直線公式
以下的公式是我從物理奧林匹亞培訓講義上抄來的。假設操縱變因（自變數）為$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}$$
<br />
## 使用 LibreOffice Calc
### 計算最接近直線斜率、截距及不準量
首先我們用 Calc 產生一些資料，假設 $y=ax+b$、$a=2$、$b=3$，再加上亂數使資料有一點偏差，產生的資料表如下：
<div style="text-align:center">
| 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 |
</div><br />
如果不想要手動輸入資料，可以直接將網頁上的資料選取並複製，直接貼到試算表中。接著我們先計算一些需要用到的數值，例如$x^2$、$y^2$、$xy$、$(\Delta x)^2$、$(\Delta y)^2$，作法相當類似，以$x^2$為例，先用滑鼠左鍵點選**儲存格F2**，接著輸入
```clike
=B2^2
```
再將滑鼠游標移到**儲存格F2**的右下角，當滑鼠游標變為黑色十字時按住滑鼠左鍵向下拖曳到**儲存格F11**，Calc會自動以公式**B3\^2**、**B4\^2**……這樣的公式計算每個儲存格對應的數值。
接著於**儲存格A14**到**儲存格J14**計算對應的數量及總合。$n$是資料點數量，先用滑鼠左鍵點選**儲存格A14**，接著輸入
```clike
=COUNT(A2:A11)
```
$\sum x$是$x$數值總合，先用滑鼠左鍵點選**儲存格B14**，接著輸入
```clike
=SUM(B2:B11)
```
以上兩個函式的詳細用法請參考官方說明書：[COUNT](https://help.libreoffice.org/Calc/Statistical_Functions_Part_One/zh-TW#COUNT)、[SUM](https://help.libreoffice.org/Calc/Mathematical_Functions#SUM)。
<br />
<img height="100%" width="100%" src="https://imgur.com/RRhGSk2.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">計算需要使用的數值</div><br />
我們可以用 Calc 內建的函式計算最接近直線的斜率，先用滑鼠左鍵點選**儲存格M2**，接著輸入
```clike
=SLOPE(D2:D11,B2:B11)
```
語法為
```clike
=SLOPE(y軸資料, x軸資料)
```
如果要手動輸入公式，先用滑鼠左鍵點選**儲存格N2**，接著輸入
```clike
=(B14*D14-A14*H14)/(B14^2-A14*F14)
```
兩者計算結果均為1.98787878787879。
<br />
再來計算最接近直線的截距，先用滑鼠左鍵點選**儲存格M3**，接著輸入
```clike
=INTERCEPT(D2:D11,B2:B11)
```
語法為
```clike
=INTERCEPT(y軸資料, x軸資料)
```
如果要手動輸入公式，先用滑鼠左鍵點選**儲存格N3**，接著輸入
```clike
=(B14*H14-D14*F14)/(B14^2-A14*F14)
```
兩者計算結果約為3.04666666666668，直到小數點後第14位才有差異。
<br />
最後計算最接近直線的$R^2$值，先用滑鼠左鍵點選**儲存格M9**，接著輸入
```clike
=RSQ(D2:D11,B2:B11)
```
語法為
```clike
=RSQ(y軸資料, x軸資料)
```
如果要手動輸入公式，先用滑鼠左鍵點選**儲存格N9**，接著輸入
```clike
=((B14*D14-A14*H14)/SQRT((B14^2-A14*F14)*(D14^2-A14*G14)))^2
```
兩者計算結果約為0.999497575579201，直到小數點後第15位才有差異。以上三個函式的詳細用法請參考[官方說明書](https://help.libreoffice.org/Chart/Trend_Lines#The_linear_regression_equation)。
<br />
<img height="60%" width="60%" src="https://imgur.com/mcs0EJs.png
" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center"></div><br />
但是 Cacl 沒有內建計算斜率、截距不準量的函數，我們只能手動輸入公式。先計算$\sigma_x$、$\sigma_y$、$\sigma$，公式分別為
```clike
=SQRT(I14/A14)
=SQRT(J14/A14)
=SQRT(N5^2+N2^2*N4^2)
```
斜率不準量$\Delta a$公式為
```clike
=SQRT(A14*N6^2/(A14*F14-B14^2))
```
計算結果為0.082813521943926。截距不準量$\Delta b$公式為
```clike
=SQRT(N6^2*F14/(A14*F14-B14^2))
```
計算結果為0.513844390399612。
<br />
### 繪製帶有誤差槓的XY散佈圖
雖然現在正式的名稱已經改為**不準量** (uncertainty) ，不過一般而言還是會稱為**誤差** (error)，在繪製XY散佈圖時資料點要加上**誤差槓** (error bar)。
#### 插入XY散佈圖
首先選取**儲存格D1:D11**（下圖1）、**B1:B11**（下圖2），再點右上角的**插入圖表**（下圖3），或是選取資料後點選單中的**插入**⇒**圖表**，插入XY散佈圖。
<img height="100%" width="100%" src="https://imgur.com/tTcmGwi.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">由工具列插入圖表</div><br />
<img height="40%" width="40%" src="https://imgur.com/ScYu57d.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">由選單插入圖表</div><br />
接著會跳出**圖表精靈**，選擇圖表類型中的**XY(散佈)**，種類為**僅限點**，再按**下一步**。
<img height="100%" width="100%" src="https://imgur.com/MVgOTD0.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">圖表精靈圖表類型</div><br />
由於我們已經先選取了資料範圍，不需要再指定資料範圍。選擇**以欄表示的為資料序列**，勾選**以第一列作為標籤**，再按**下一步**。
<img height="100%" width="100%" src="https://imgur.com/tjM4Gi0.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">圖表精靈資料範圍</div><br />
由於我們目前只想在這張圖中畫出一組資料序列，直接按**下一步**即可。
<img height="100%" width="100%" src="https://imgur.com/WpoxVz1.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">圖表精靈資料序列</div><br />
我們習慣上不會加上**題名**和**副題**，會標示在圖的說明文字中，這兩格空白即可。**X軸**、**Y軸**這兩格寫上代表的物理量及單位。由於圖上只有一組數據，所以不需要顯示圖例。格線的部分則是依照個人喜好決定，沒有一定的作法。按**完成**之後產生的圖形如下。
<img height="100%" width="100%" src="https://imgur.com/ggAeg9D.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">圖表精靈圖表元素</div><br />
<img height="70%" width="70%" src="https://imgur.com/w8eqHyt.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">只有資料點的XY散佈圖</div><br />
#### 修改圖片格式
接下來還可以再修改資料點的格式，先在圖上連點按滑鼠左鍵兩下進入圖表，再於資料點上按滑鼠右鍵開啟快速選單，選取**設定資料序列格式**。
<img height="70%" width="70%" src="https://imgur.com/K1fvSDL.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">快速選單</div><br />
在資料序列視窗中的**線條**分頁，可以修改線條的樣式及顏色，也可以修改資料點圖示。
<img height="70%" width="70%" src="https://imgur.com/iShD4m8.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">資料序列視窗</div><br />
從右側的**選取**⇒**圖庫**，選取自己喜歡的樣式。
<img height="100%" width="100%" src="https://imgur.com/EOvcf4Y.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">資料序列圖示</div><br />
<img height="70%" width="70%" src="https://imgur.com/uQBPjsY.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改資料點圖示的XY散佈圖</div><br />
如果覺得預設的座標軸標示、標籤的字體太小，可以在軸上面按滑鼠右鍵，從選單中選取**設定軸格式**。
<img height="70%" width="70%" src="https://imgur.com/rMbRKnu.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">設定軸格式</div><br />
選擇自己喜歡的字型、樣式及大小，再按下**確定**即可。
<img height="80%" width="80%" src="https://imgur.com/AMSLiyT.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">設定Y軸字型</div><br />
我們可以用類似的方法修改座標軸標籤。
<img height="70%" width="70%" src="https://imgur.com/Q4gqmMw.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改座標軸標籤</div><br />
<img height="80%" width="80%" src="https://imgur.com/vYfiVSg.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">設定Y軸標題字型</div><br />
<img height="70%" width="70%" src="https://imgur.com/gMm7p21.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改座標軸格式和標籤後的XY散佈圖</div><br />
#### 插入誤差槓
在資料點上按滑鼠右鍵，從選單中選取**插入X誤差線**。
<img height="70%" width="70%" src="https://imgur.com/EBKXmYY.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">插入X誤差線選單</div><br />
由於我們的C欄資料就是$x$軸不準量，所以我們要選擇**誤差類別**中的**儲存格範圍**（下圖1），再按下圖2的標籤從資料表中選取儲存格C2:C11，由於我們的不準量正負值相同，只要勾選下圖3的**兩者同值**，最後按下**確定**（下圖4）即可。接著可以用類似的方法插入Y誤差線。
<img height="70%" width="70%" src="https://imgur.com/jc9iKhd.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">資料序列「y」的X誤差線</div><br />
<img height="70%" width="70%" src="https://imgur.com/5nPv5Zc.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">插入誤差線後的XY散佈圖</div><br />
#### 插入趨勢線
如果要在圖上加入趨勢線，可以在資料點上按滑鼠右鍵，從選單中選取**插入趨勢線**。在下圖的視窗中選取**線性**（下圖1），如果要在圖上**顯示方程式**及**決定係數**$R^2$，則要勾選圖中標籤2、3，最後再按下**確定**（下圖4）即可。但是我們通常不會在圖上顯示方程式及決定係數，會另外標示於圖表說明中。
<img height="50%" width="50%" src="https://imgur.com/3f0aqKG.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">資料序列「y」的趨勢線</div><br />
<img height="70%" width="70%" src="https://imgur.com/yOJHyBC.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">插入擬合線後的XY散佈圖</div><br />
#### 匯出成圖檔
如果要把XY散佈圖匯出成圖檔，不要進入編輯圖表的模式，直接於圖上按滑鼠右鍵，從選單中選取**匯出為影像**，決定存檔路徑、檔名及圖片格式，最後按**存檔**。
<img height="80%" width="80%" src="https://imgur.com/QKLfofC.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">匯出為影像選單</div><br />
<img height="80%" width="80%" src="https://imgur.com/J5rSMdp.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">決定存檔路徑、檔名及圖片格式</div><br />
<img height="70%" width="70%" src="https://imgur.com/ZxVFm0M.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">匯出後的影像</div><br />
## 使用 SciDAVis
### 從 CSV 檔匯入資料
先將以下的資料複製、貼上到文字編輯器中，將檔案儲存為**data.csv**，文字編碼應該用 ASCII 或 UTF-8 皆可。
```python
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
```
<br />
由選單：**File**⇒**Import ASCII** 或 **Ctrl + K**，開啟匯入檔案視窗。
<img style="display: block; margin-left: auto; margin-right: auto" height="80%" width="80%" src="https://imgur.com/wpyHIFj.png">
<div style="text-align:center">匯入檔案選單</div><br />
由於資料是以逗號分隔，在匯入檔案視窗中，要記得**將 Separator 選項改成逗號**。另外也可以選擇要將資料匯入哪一個表格，可以蓋掉現有表格的資料，也可以另外增表格。
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/iCZeuWg.png">
<div style="text-align:center">匯入檔案視窗</div>
<br />
<img style="display: block; margin-left: auto; margin-right: auto" height="100%" width="100%" src="https://imgur.com/udvWgbU.png">
<div style="text-align:center">匯入資料檔後的表格</div>
<br />
### 設定欄位資料格式
接著在要修改的欄位上按滑鼠左鍵，在畫面右側也可以更改欄位的描述及資料格式。
<img height="40%" width="40%" src="https://imgur.com/LohW49w.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">更改欄位的描述</div>
<br />
<img height="40%" width="40%" src="https://imgur.com/mTs4Em1.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">更改欄位的資料格式</div>
<br />
在要修改的欄位上按滑鼠右鍵，選擇**Set Column(s) As**即可修改欄位的資料用途，分別將欄位index設為X、欄位x設為X、欄位errx設為X Error、欄位y設為Y、欄位erry設為Y Error。
<img height="80%" width="80%" src="https://imgur.com/7HUNDw9.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">設定欄位的資料用途</div>
<br />
### 繪製帶有誤差槓的XY散佈圖
#### 插入XY散佈圖
先選取欄位x、y之後，再從選單選取**Plot**⇒**Scatter**，或是在欄位上按滑鼠右鍵，從快速選單中選取**Plot**⇒**Scatter**。
<img height="90%" width="90%" src="https://imgur.com/20DdYJZ.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">從選單插入XY散佈圖</div>
<br />
<img height="80%" width="80%" src="https://imgur.com/9MtY8Ph.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">從快速選單插入XY散佈圖</div>
<br />
<img height="60%" width="60%" src="https://imgur.com/O0ZKIXS.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">預設格式的XY散佈圖</div>
<br />
#### 修改圖片格式
如果想要修改資料點格式，可以在圖片的資料點上連點滑鼠左鍵兩下，勾選**Fill Color**會在點中填滿指定的顏色，**Edge Color**是點的邊緣線條顏色。
<img height="90%" width="90%" src="https://imgur.com/R8XlK61.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改資料點格式</div>
<br />
點選**Plot Associations**，可以修改資料點x軸、y軸數值對應的欄位。
<img height="40%" width="40%" src="https://imgur.com/t1sgbnz.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改資料點對應的欄</div><br />
如果想要修改座標軸及標籤的格式，可以在圖片的座標軸上連點滑鼠左鍵兩下，點選**Axis**分頁，在**Title**欄位中可以修改標籤內容，支援HTML語法；點選**Axis Font**可以修改座標軸數字的字型及大小。
<img height="80%" width="80%" src="https://imgur.com/Hr27J0S.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改座標軸性質</div>
<br />
<img height="60%" width="60%" src="https://imgur.com/aUHG1m0.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改座標軸數字的字型</div>
<br />
<img height="60%" width="60%" src="https://imgur.com/uThGO80.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">修改格式後的XY散佈圖</div>
<br />
#### 插入誤差槓
點選圖片後從選單選取**Graph**⇒**Add Error Bars**或按快速鍵**Ctrl+B**。
<img height="40%" width="40%" src="https://imgur.com/0z3DTnS.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">插入誤差槓</div>
<br />
勾選**Existing column**，從資料表中選取不準量的來源，在視窗的下方勾選**X Error Bars**或**Y Error Bars**，最後點選右上方的**Add**。
<img height="60%" width="60%" src="https://imgur.com/iJxgKXg.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center"></div>
<br />
<img height="60%" width="60%" src="https://imgur.com/8qUo9yn.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">插入誤差槓的XY散佈圖</div>
<br />
#### 插入趨勢線
點選圖片後從選單選取**Analysis**⇒**Quick Fit**⇒**Fit Linear**。
<img height="55%" width="55%" src="https://imgur.com/dI4rGZ3.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">選單插入趨勢線</div>
<br />
擬合的結果如下，看起來只有考慮y軸不準量。
$$ a = 1.9979669965615 \pm 0.0715310360117212 $$
$$ b = 2.99177166800124 \pm 0.371842523931618 $$
$$ R^2 = 0.999509982857237 $$
<img height="70%" width="70%" src="https://imgur.com/RRg1jNQ.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">擬合結果</div>
<br />
如果要再開啟輸出結果的視窗，可以從選單選取**View**⇒**Results Log**。
<img height="50%" width="50%" src="https://imgur.com/gLZL2OS.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">開啟輸出結果視窗</div>
<br />
#### 匯出成圖檔
點選圖片後從選單選取**File**⇒**Export Graph**⇒**Current**或按快速鍵**Alt+G**。
<img height="50%" width="50%" src="https://imgur.com/UGbQt8A.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">匯出圖檔選單</div>
<br />
決定存檔路徑、檔名及圖片格式，最後按**存檔**。經過測試後發現，存成png檔大約400KB，存成jpg檔大約35KB，但是在電腦螢幕上看起來效果差不多。
<img height="70%" width="70%" src="https://imgur.com/Dlc2N02.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">匯出圖檔視窗</div>
<br />
<img height="60%" width="60%" src="https://imgur.com/f8m3sYL.jpg" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">匯出的XY散佈圖</div>
<br />
## 使用 Python
最後我們用Python處理數據，需要用到的函式庫為**NumPy**、**SciPy**及**Matplotlib**。程式碼如下
```python=
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()
```
<br />
### 從 CSV 檔匯入資料
使用的函式為**numpy.loadtxt**，詳細的用法請參考[官方說明書](https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html)。程式碼第6行
```python
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。
<br />
### 方法1: 手動輸入公式
我們先手動算出需要用到的數值。例如程式碼第15行，用函式**len**計算資料數量再指定物件n。語法有兩種
```python
[陣列物件].len
len([陣列物件])
```
程式碼第16到22行，用函式**numpy.sum**計算陣列資料的總合，由於資料為1維陣列語法為
```python
[陣列物件].sum
np.sum([陣列物件])
```
詳細的用法請參考[官方說明書](https://docs.scipy.org/doc/numpy/reference/generated/numpy.sum.html)。
<br />
程式碼第26到30行，將數值代入文章最前面的公式，計算結果為
$$ a = 1.987878787879 \pm 0.082813521944 $$
$$ b = 3.046666666667 \pm 0.513844390400 $$
$$ R^2 = 0.999497575579 $$
<br />
### 方法2: 使用內建函數
我們使用**scipy.optimize.curve_fit**計算最接近直線的斜率、截距及不準量，由於程式一開始已經引入這個函式，所以語法為
```python
curve_fit([擬合函數型式], [x軸資料], [y軸資料], sigma=[y軸不準量], absolute_sigma=True)
```
我們在程式中寫成
```python
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軸資料的不準量。詳細的用法請參考[官方說明書](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)。
<br />
接著使用**numpy.corrcoef**計算相關係數$R$，語法為
```python
np.corrcoef([x軸資料], [y軸資料])
```
回傳值是一個2*2的陣列，我們要的$R^2$陣列右上角元素值的平方，因此程式碼寫為
```python
cor = np.corrcoef(x, y)[0, 1]
print("r-squared: %.12f" % cor**2)
```
詳細的用法請參考[官方說明書](https://docs.scipy.org/doc/numpy/reference/generated/numpy.corrcoef.html)。
<br />
計算結果為
$$ a = 1.997966995944 \pm 0.071531035658 $$
$$ b = 2.991771672485 \pm 0.371842522314 $$
$$ R^2 = 0.999497575579 $$
<br />
### 繪製帶有誤差槓的XY散佈圖
最後我們用matplotlib.pyplot繪圖，由於函式庫的全名太長，通常引入時會寫成
```python
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，詳細的用法請參考[官方說明書](https://matplotlib.org/3.1.1/api/_as_gen/matplotlib.pyplot.errorbar.html)。
7. 第65、66行：為了加上最接近直線，我們先用擬合得到的斜率**popt[0]**、截距**popt[1]**，將陣列x資料代入計算對應的y值，再用plot函式畫線。
<img height="60%" width="60%" src="https://imgur.com/blWlHKn.png" style="display: block; margin-left: auto; margin-right: auto;"/>
<div style="text-align:center">用 matplotlib 繪製的XY散佈圖</div><br />
## 結語
以下是不同工具的計算結果，其中SciDAVis和Python(內建)在擬合時會考慮y軸不準量，而且計算結果相當接近，如果考慮學生學習工具上的難度，我會建議學生使用SciDAVis，如果不排斥看到程式碼再學著使用Python。我相信按照這個標準得到的XY散佈圖，應該可以符合科展評審的標準。
<div style="text-align:center">
| | 內建公式 | 手動輸入 | 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 |
</div>
---
###### tags:`Physics`、`LibreOffice`、`SciDAVis`、`Python`、`校訂必修`
Sign in with Wallet
Connect another wallet