---
# System prepended metadata
title: 【R 語言與統計資料分析】Ch9：推論統計與中央極限定理
tags: [R 語言與統計資料分析]
---
---
title: 【R 語言與統計資料分析】Ch9：推論統計與中央極限定理
image: https://ppt.cc/fnT4xx@.png
---
# 【R 語言與統計資料分析】Ch9：推論統計與中央極限定理
樣本與母體間的關係
讓我們回到數據分析的核心概念，我們希望透過已知的部分資訊 (樣本) 去「猜測」真實的母體狀態，而這個猜測的方法就稱為「統計推論」。
# 9.1 從估計開始
若我們今天收到一筆樣本，是學生在問卷上填寫的身高資料，而我們想透過此資料進行健康政策研究，那們要如何才能做出是當的政策判斷依據呢? 以下是美國學生的身高數據，那我們先換成公分，比較符合台灣人的直覺。
```r
library(dslabs) #載入套件中的資料集
# 如果沒有此套件，請執行 install.packages("dslabs")
data(heights)
# 觀看查詢前 6 筆資料
head(heights)
```
heights 資料集內有性別與身高 (英吋) 資料
```r
heights$height_cm <- heights$height * 2.54 #英吋換成公分
# 觀看查詢新版前 6 筆資料
head(heights)
```
heights 資料集內有性別與身高 (公分) 資料
接著，讓我們繪製此資料集的直方圖，也就是實證機率密度函數 (empirical probability density function)，來看看資料的分佈狀況。
```r
hist(heights$height_cm)
```
height (cm) 的直方圖
我們可以發現，身高確實符合鐘形分配，但接下來還會面對許多問題:
- 我們真的可以用這組樣本代表母體嗎?
- 這個母體真的是常態分佈嗎?
- ... ...
在這個章節，我們會用統計學的思維以及 R 語言的實證回答這些問題。
## 9.1.1 點估計
凡事總是要有個起點，我們首先從單點估計開始說起，若想了解一個資料的分布型態，我們當然會想要了解其他關鍵指標，比如說平均數與標準差，因為當我們找到這些指標後，就可以對母體有更清楚的描述，而針對單一指標的估計，就被稱為「**點估計**」(point estimator)。
::: success
💡 **樣本動差**
在統計學中，「動差」（Moments）是用來描述分佈形狀的數學特徵，我們在此介紹幾個重要動差:
- **平均數** (一階動差): 為對母體期望值的估計。
$$\bar{x} = \frac{1}{n} \sum_{i=1}^{n} x_i$$
- **變異數** (二階中央動差): 描述數據的離散程度，在大樣本的情況之下，分母為 $n$ 或 $n-1$ 差異不大。
$$s^2 = \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2$$
- **偏態** (Skewness，三階標準化動差): 描述分佈的左右對稱性。
$$g_1 = \frac{\frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^3}{\left[ \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 \right]^{3/2}}$$
- $g_1 > 0$：右偏（右尾較長，資料集中左側）。
- $g_1 < 0$：左偏（左尾較長，資料集中右側）。
- **峰態** (Kurtosis，四階標準化動差): 描述分佈頂部的尖銳度及尾部的厚度。
$$g_2 = \frac{\frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^4}{\left[ \frac{1}{n} \sum_{i=1}^{n} (x_i - \bar{x})^2 \right]^2}$$
- $g_2 = 3$：常態分配（基準值）。
- $g_2 > 3$：高狹峰（Leptokurtic），尾巴比常態分配厚（極端值多）。
- $g_2 < 3$：低闊峰（Platykurtic），分佈較平緩。
接著，讓我們用 R 來計算這些指標:
```
#平均數
mean(heights$height_cm)
#變異數
var(heights$height_cm)
# 安裝與載入套件
# install.packages("moments")
library(moments)
#偏態
skewness(heights$height_cm) # 若接近 0 則代表對稱
#峰態
kurtosis(heights$height_cm) # 若接近 3 則代表符合常態分配
```
:::
## 9.1.2 區間估計
在了解點估計後，接下來要回答的問題就是: 點估計準嗎? 一個簡單的答案是: 「**肯定不準**」，是不太可能直接猜中母體動差的，但是我們是可以用一個區間來衡量猜測的準確度，而這就是**區間估計** (interval estimation)!
以母體均數為例，我們希望找到一條式子，滿足以下條件，而此處的上下界所包含的範圍，就是區間估計。
$$P[L \leq \mu \leq U]=95\%$$
在理論上，我們必須先建立一個包含要估計的參數的樞軸量 (pivotal)，但其抽樣分配又與未知參數無關，且我們必須知道它的抽樣分配 (不然無從算起)，我們在此處就先跳過不談細節，若對其數學方法有興趣，可以參考以下連結之[網路文件](https://en.wikipedia.org/wiki/Pivotal_quantity)。
以常態分配為例，我們所計算之 95% 信賴區間為 (假設母體標準差已知):
$$\left[ \bar{X} - 1.96 \frac{\sigma}{\sqrt{n}}, \quad \bar{X} + 1.96 \frac{\sigma}{\sqrt{n}} \right]$$
若我們算出 $\mu$ 的 95% 信賴區間為 $[10,20]$，是代表若我們重複此隨機試驗 100 次，會約有 95 次所算出來的區間會確實涵蓋到母體均數，而非此次抽樣的結果涵蓋真正均數的機率為 95%，因為現實只會有「**沒涵蓋到**」與「**有涵蓋到**」兩種狀況，因此機率不是 1 就是 0。
信賴區間的含意
圖片來源:林惠玲、陳正倉(2004) 統計學:方法與應用
# 9.2 弱大數法則與中央極限定理
## 9.2.1 弱大數法則
弱大數法則描述隨樣本數 $n$ 增加，樣本估計與母體參數的差距而趨近於 0，簡單來說就是當樣本數變多時，資料會越準。
$\lim _{n\rightarrow \infty }P(|\bar{X}_{n}-\mu |>\epsilon )=0$
```r
# 探討樣本數與均數間的關係
set.seed(888888)
n <- c(10,30,50,100,1000,10000,100000) #抽樣次數
average <- c() #建立空白向量
for(i in 1:length(n)){
average[i] <- mean(rpois(n[i],50)) #卜瓦松分配 lambda = 50 也就是平均數
}
df <- data.frame(n = n,average = average, diff=(50-average))
print(df)
```
## 9.2.1 中央極限定理
在現實世界中，我們基本上不太可能知道資料真實的分佈情形，那此時就無法像前述一樣簡單地計算出樞軸量 (pivotal)。但慶幸的是，當樣本數趨近無窮大時，**樣本平均數之抽樣分佈會收斂**(convergence in distribution)至**常態分佈**，因此對於母體均數而言，我們還是可以做出適當的統計推論的！
$$\frac{\bar{X}_n-\mu}{\frac{\sigma}{\sqrt{n}}}\xrightarrow[n\rightarrow\infty]{d}N(0,1)$$
(我們可用樣本標準差 $s$ 來取代 $\sigma$)
#### 以卡方分配為例
```r
n = c(1, 10, 100)
color = c("red", "blue", "green")
nSim <- 10000
for(i in 1:length(n)){
# 一次生成矩陣：列數為模擬次數，行數為樣本數
# 然後計算每一列的平均值，結果就是一個向量
meanVector <- rowMeans(matrix(rchisq(nSim * n[i], df=5), nrow = nSim))
plot(density(meanVector), col = color[i], xlim=c(0, 15),
ylim=c(0, 1.3), main = "CLT on chisquare dist.(df=5)", xlab = "x-bar")
par(new = TRUE)
}
legend(x="topright", legend=c("n=1", "n=10", "n=100"),
lty = 1, col = color, lwd = 1.5)
```
卡方分配的樣本平均收斂速度
#### 以二項分配為例
```
n = c(1, 30, 100)
color = c("red", "blue", "green")
nSim <- 10000
p_success <- 0.2 # 設定二項分配的成功機率
for(i in 1:length(n)){
# rbinom(n, size, prob):
# 這裡 size=1 代表白努力試驗，我們抽 nSim * n[i] 個，再排成矩陣算平均
raw_data <- rbinom(nSim * n[i], size = 1, prob = p_success)
meanVector <- rowMeans(matrix(raw_data, nrow = nSim))
# 繪圖設定：xlim 改為 0 到 1 之間，因為二項分配平均值（機率）範圍在此
plot(density(meanVector), col = color[i], xlim=c(0, 1),
ylim=c(0, 15), main = "CLT on Binomial dist. (p=0.2)", xlab = "x-bar (Sample Proportion)")
par(new = TRUE)
}
legend(x="topright", legend=c("n=1", "n=30", "n=100"),
lty = 1, col = color, lwd = 1.5)
```
二項分配的樣本平均收斂速度較慢
#### 以對數常態分配為例
對數常態分配，為隨機變數 $X$ 的「自然對數值」$\ln(X)$ 服從常態分配。也就是，$\ln(X) \sim N(\mu, \sigma^2)$。
```
n_list = c(1, 10, 100, 500, 10000)
color = c("red", "blue", "green", "purple","orange")
nSim <- 10000
for(i in 1:length(n_list)){
# 對數常態分配 rlnorm
meanVector <- rowMeans(matrix(rlnorm(nSim * n_list[i], meanlog=0, sdlog=2), nrow = nSim))
plot(density(meanVector), col = color[i],
main = "CLT on Log-normal (Very Slow Convergence)", xlab = 'xbar',
xlim = c(0, 50), ylim = c(0, 0.5))
par(new = TRUE)
}
legend("topright", legend=paste("n=", n_list), col=color, lty=1)
```
對數常態分配的樣本平均收斂速度極緩慢且會出現右長尾
\
由這些案例可知，中央極限定理雖然理論上成立，但是根據不同的分配會有不同的收斂速度。因此，我們需要對資料有更深的認識才能做出正確的統計推論。 (這也是為何我們總是要先繪製資料的直方圖)
# 9.3 假設檢定
我們在做統計推論時，重點就是為了「回答問題」，那回答問題的一種方式就是**假設檢定**。
透過假設檢定，我們可以試著回答這些問題
- 建中男同學都不高，只有 173 公分，真的假的?
- 哪一個手機軟體介面會讓用戶有更好的付費意願?
- 在 50 公尺自由式比賽中，被分到 5~8 水道 的選手成績普遍比 1~4 水道的選手好 ?
## 9.3.1 虛無假設與對立假設
在假設檢定中，我們會建立兩個互斥的假設，分別為：
- **虛無假設** (null hypothesis)：這是我們要檢定的假設，一般以 $H_0$ 表示。
簡單來說，我想確認母體平均身高是否可能為 173 公分，那此時的 $H_0$ 就設定為： $H_0 : \mu = 173$
- **對立假設** (alternative hypothesis)：這是與虛無假設相悖的假設，以 $H_1$ 或 $H_A$ 表示。
承上例：$H_1： \mu > 173$
若透過檢定認為**虛無假設不成立**，則我們會反過頭來**接受對立假設**。
透過檢定，我們會獲得下列兩個結果：
- 拒絕 $H_0$ 並接受 $H_1$。
- 無法拒絕 $H_0$。
::: info
💡 **假設檢定的基本想法**
若我們先**假定虛無假設為真**，我們抽取出來的樣本就不應該偏離這個值太多，**若出現極端值**的情形，則我們有**信心說我們拒絕虛無假設**，而若沒有出現這種極端的情形，則我們只能說我們**無法拒絕虛無假設**。
:::
## 9.3.2 顯著水準與拒絕域
#### 顯著水準
那麼出現這種極端值的機率是多少呢？一般來說，我們會定義**顯著水準** (significance level) 為 $\alpha$ ，並令 $\alpha=0.01, 0.05\ or\ 0.1$。顯著水準也就是**當** $H_0$ **成立時**，**發生某極端事件的機率**。不同的顯著水準也會影響檢定被拒絕與否。
而要計算機率，那我們就得建構出我們有興趣的統計量的「機率分配」，此部分就會與「**虛無假設**」有密切關係，因為我們進行假設檢定時，會先**假定虛無假設為真**，此時所算出來的機率就是使用**虛無分配** (null distribution) 去做計算。
#### 拒絕域
所謂拒絕域 (rejection region, RR) 是在虛無假設成立之下，並給定機率 $\alpha$，所圍出來的範圍，而這個範圍會有一個臨界值 (critical value)，當我們**抽出來的估計實現值**超過臨界值而落在拒絕域時，我們就可以有信心地去拒絕虛無假設。
>注意：臨界值是經由虛無分配所堆導出之值，相對於抽出的樣本資料 (ex post)，是一個「事前 (ex ante)」的概念
我們可以回到身高的例子，$H_0 : \mu = 173$ 且 $H_1： \mu > 173$。那抽出來的樣本平均到底要要超過多少，我們才會拒絕虛無假設呢？是 175、180、還是 200 公分呢？
所以我們的目的就是找出一個臨界值 $c$ 使得下式成立：
$P(\bar{X} \geq c |\mu=173)=\alpha$
也就是說，當我們所抽取出的**樣本平均實現值 $\bar{x} \geq c$** 時，就是在 $H_0$ 為真的情況之下，真的出現了極端值，那我們就可以**在顯著水準 $\alpha$ 之下拒絕虛無假設**。
虛無分配的實際呈現
## 9.3.3 假設檢定的步驟
在假設檢定中，我們對**母體**分配的「**參數**」是有興趣的，如常態分配中的 $\mu$ ，或是 Bernoulli 分配中的 $p$。首先，我們其實不知道母體的真實分配，所獲得的資料又是抽出的**樣本**，那這樣子獲得的樣本資料有辦法去推論母體嗎？我們可以透過一些**統計理論**，如中央極限定理等等，協助我們進行參數的檢定。接著，我們來介紹假設檢定的基本步驟，我們以 $\theta$ 來表示想要估計的參數， $\theta^*$ 表示我們在**虛無假設**下該參數的值。
1. 設定虛無假設 $H_0$ 以及對立假設 $H_1$
1. $H_0: \theta= \theta^*$
2. $H_1$ 根據我們所獲得資訊的不同或是有特別的檢定目的，可以分為**雙尾檢定**以及**單尾檢定：**
1. 雙尾檢定(Two-Tailed Test)：$\theta \neq \theta^*$
2. 右尾檢定(Right-Tailed Test)：$\theta > \theta^*$
3. 左尾檢定(Left-Tailed Test)：$\theta < \theta^*$
圖片擷取自：D.Bang，用人話說明什麼是假設檢定
2. 建立樞紐統計量(pivotal quantity) $\varphi(\theta,X_i)$，我們透過找出 $\varphi(\bullet)$的抽樣分配來進行統計推論。
:::success
$\varphi(\theta,X_i)$ 中沒有任何未知參數，並且它的分配與 $\theta$ 無關。
若$\varphi(\bullet)$的抽樣分配未知，我們可以透過**中央極限定理**以 $N(0,1)$ 來近似。
中央極限定理：$\frac{\bar{X_n}-\mu}{\sqrt{Var(\bar{X_n})}}\xrightarrow[]{d}N(0,1)$
:::
3. 選擇顯著水準 $\alpha$ ，傳統上我們會設定 $\alpha=0.01,0.05,0.1$。
4. 透過$\varphi(\bullet)$來找出對應顯著水準 $\alpha$ 的**臨界值 $\varphi^*$**，並以此建構**拒絕域**。(不同檢定各自的拒絕域的位置可以對照上圖分配之著色處)
5. 檢視實際抽出的資料 $\varphi_0=\varphi(x_1,x_2,...,x_n,\theta_0)$，是否落在**拒絕域**中。當我們的**統計量掉在拒絕域**中的話，我們就可以**拒絕虛無假設**。
#### 假設檢定與信賴區間
假設檢定與信賴區間為一體之兩面，在 null distribution 成立之下，我們也可以建立包含母體參數之 $(1-\alpha)\%$信賴區間，如果檢定統計量落在這個區間之外，則可拒絕 $H_0$。
:::info
📚 **民意調查的案例**
我們隨機調查 1000 名受訪者，詢問他們喜不喜歡吃香菜，是的話紀錄為$1$，否的話就紀錄為$0$。因此可將本題受訪者偏好的樣本視為從 $Bernoulli(p)$ 中抽出的 $i.i.d$ 隨機變數。
$i.e$ $\left \{ X_i \right \}_{i=1}^{1000} \sim^{iid}Bernoulli(p)$，其中
$X_i={\begin{cases}1,\ 愛吃香菜 \\0,\ 不愛吃香菜 \end{cases}},\ i=1,2,...,1000$
我們發現有 570 位受訪者覺得香菜很讚，有點介在中間啊，故我們想要知道 $p$ 是否為 $0.5$，我們設定顯著水準為 $0.05$。
:::
Solution:
1. 我們首先建立 $H_0$ 與 $H_1$。
$H_0: p = 0.5$
$H_1: p \neq 0.5$
2. 建立樞紐量 pivot，我們預估計的參數是而我們是用 $\bar{X}$ 去估計 $p$，但由於我們不知道 $\bar{X}$ 的分配，故我們透過**中央極限定理**去幫助我們建立樞紐量。
$\frac{\bar{X_n}-E(\bar{X_n})}{\sqrt{Var(\bar{X_n})}} = \frac{\bar{X_n}- p}{\sqrt{\frac{p(1-p)}{n}}}\xrightarrow[]{d}N(0,1)$
3. 建立 $\alpha=0.05$ 的**拒絕域**。由於 $N(0,1)$ 是對稱的，在雙尾檢定之下，我們只要找到 $1-\frac{\alpha}{2}$ 的 quantile 即可，我們可以透過 R 語言，也可以透過傳統的查表來找到這個數值
在 R 中，我們可使用寫好的套件直接進行假設檢定:
```r
# 設定數據
n <- 1000 # 樣本總數
x <- 570 # 覺得很讚的人數
p_null <- 0.5 # 虛無假設的比例
# 方法一：精確二項檢定 (Exact Binomial Test)
res_exact <- binom.test(x, n, p = p_null, alternative = "two.sided")
print(res_exact)
# 方法二：比例檢定 (Proportion Test, 使用常態近似)
res_prop <- prop.test(x, n, p = p_null, alternative = "two.sided", correct = TRUE)
print(res_prop)
```
Exact binomial test 檢定報表
## 9.3.4 檢定的 p 值
在 $H_0$ 假設為真的情形之下，比我們所觀察到之觀測值更極端之機率。
因此，若 p-value < $\alpha$，我們就可以拒絕 $H_0$。也就是在虛無假設成立的情況之下，我們竟然觀察到這麼低機率的事件，因為太誇張了，所以我們拒絕無假設。
$\bullet$ 拒絕 $H_0$：p-value $\leq \alpha$
$\bullet$ 無法拒絕 $H_0$：p-value $\geq \alpha$
讓我們回到前述民調的例子，我們可以透過下列式子來計算 p-value：
$p\ value=P(|\bar{X_n}| \geq 0.57|p=0.5)=P(|\frac{\bar{X_n}- 0.5}{\sqrt{\frac{0.5(1-0.5)}{1000}}}|
\geq \frac{0.57 - 0.5}{\sqrt{\frac{0.5(1-0.5)}{1000}}})=\\P(|Z| \geq 4.43)=2 \times P(Z \geq 4.43)=9.423 \times 10^{-6}$
(注：後面有 $\times 2$ 是來自於雙尾檢定拆絕對值，若是單尾檢定就不用 $\times 2$ )
我們可以看出 p-value 超級小，因此我們也就**有信心拒絕虛無假設**。
## 9.3.5 檢定也會犯錯 - 型一誤差與型二誤差
我們在搜集資料，然後進行統計資料，是否真有可能就那麼衰，真的抽樣抽到了極端值，導致我們錯誤拒絕了 $H_0$ 呢？是否也有可能錯誤地無法拒絕 $H_0$ ，但其實 $H_1$ 才是對的呢？
在此，我們以下表表示上面的情形：
| | 無法拒絕 $H_0$ | 拒絕 $H_0$ |
| --- | --- | --- |
| $H_0$ 是正確的 | 決策無誤 | $Type\ I\ error$ |
| $H_0$ 是錯誤的 | $Type\ II\ error$ | 決策無誤 |
- 型一誤差( $Type\ I\ error$ )：錯誤拒絕 $H_0$，錯殺無辜。
- 型二誤差( $Type\ II\ error$ )：錯誤接受(無法拒絕) $H_0$，錯放壞人。
在檢定中我們會希望**控制型一誤差與型二誤差發生的機率**，當然是希望越小越好，只可惜**魚與熊掌不可兼得**，若樣本數 $n$ 不變的情況之下，兩者存在 trade-off 的關係。
- 犯下型一誤差的機率：$\alpha = P(型一誤差)=P(拒絕H_0|H_0為真)$
這個 $\alpha$ 在簡單假設的情況之下，**就是我們的顯著水準**， $\alpha$ 也可以說是我們可容忍犯型一誤差的最大機率。
> 是否有機會「洗」p 值 ?
- 犯下型二誤差的機率：$\beta= P(型二誤差)=P(無法拒絕H_0|H_0為假)$
而因為 $H_0$ 是有問題的，我們可以知道 $\beta$ 的機率來自於母體**未知**參數的真值，在不同的 $\theta$ 的設定之下也會有不同的機率值，故 **$\beta$ 是 $\theta$ 的函數 ( $\beta=\beta(\theta)$ )。**
:::success
💡 **檢定力函數** (Power Function)
我們也因此可以延伸出一個觀念 ——「**檢定力函數**」，定義檢定力函數如下：
$\pi(m)=P(拒絕H_0|\theta=m)=1-\beta(m)$ ← 想想看，why?
也就是說，若 $H_1$ 為真，**檢定力就是不犯型二誤差的機率**。
:::
型一誤與型二誤差
# 9.4 檢定不顯著怎麼辦
檢定不顯著可能只是代表目前的數據不足以佐證我們的猜想，因此我們還無法推翻虛無假設。另一方面，也要去思考是否研究設計有問題? 或是我的對立假設真的是錯的?
## 9.4.1 如何降低 p 值
1. **增加樣本數** ($n$)
我們以中央極限定理為例，樣本數越大，標準誤（Standard Error）就越小，此時就有機會在雷同的情境下拒絕虛無假設。
$$CI = \bar{X} \pm Z_{\alpha/2} \cdot \frac{\sigma}{\sqrt{n}}$$
但要小心「過度顯著」，當 $n$ 極大時，微不足道的差異也會顯著。
2. **減少測量誤差與雜訊** (Noise Reduction)
如果你無法增加樣本，那就設法讓現有的數據更「精準」。實驗控制：在更嚴格受控的環境下進行實驗，減少環境干擾。精進工具：使用精確度更高的儀器。同樣地，如果你知道某些因素（如年齡、性別）會影響結果，將它們納入模型作為控制變數，可以移除「背景雜訊」，從而讓主效應的 P 值下降。
## 9.4.2 請不要「操弄」p 值
1. **操弄統計樣本**:
研究者在收集數據的過程中，每收幾筆就跑一次檢定：如果 $P > 0.05$，就繼續收數據。一旦 $P < 0.05$，就立刻停止實驗並宣稱結果顯著。
這會讓第一型錯誤 ($\alpha$) 遠高於預設的 0.05，因為你給了數據多次「偶然達標」的機會，雷同於我們之前說的「洗」p 值。
2. **選擇性報告** (Selective Reporting / Cherry Picking):
一個實驗可能測量了 20 種不同的指標（例如：血壓、心跳、睡眠、食慾等）。如果只有「食慾」達到了 $P < 0.05$，研究者就只報告食慾的結果，並假裝沒測量過其他 19 項。後果：根據機率，同時測量 20 項指標時，純粹靠運氣出現一項顯著結果的機率高達 $64\%$。
3. **刪除極端值** (Outlier Manipulation):
人為地定義哪些數據是「異常值」，但極端值不一定是異常值。研究者會嘗試刪除某些數據點，看看 $P$ 值會不會降下來。如果刪除後變顯著，就找個藉口（如「儀器誤差」）把該點踢除。
4. **變更分組方式** (Data Binning):
將連續變數人為地切割成類別變數。例如：嘗試將年齡切成「20-30, 31-40...」，如果沒效，就改切成「25歲以下, 25歲以上」，直到跑出顯著結果為止。
# 參考資料
1. 陳旭昇(2024)，資料分析的統計學基礎:使用R語言，東華書局
2. 林建甫 Jeff Lin(2020)，[R 資料科學與統計](https://bookdown.org/jefflinmd38/r4biost)
3. Gareth James, Daniela Witten, Trevor Hastie, Robert Tibshirani. (2017). An Introduction to Statistical Learning: With Applications in R. New York: Springer.
4. 陳基國(2024). 基礎統計與R語言. 台北：五南圖書出版股份有限公司