社會科學統計方法

統計估計、假設檢定

1 課程目標

我們用樣本統計量來估計母體參數，例如平均數、標準差。統計估計可以分成點估計以及區間估計。

本週上課將介紹統計估計的意義，如何建立區間估計，以及如何檢證連續變數的平均值是否不等於0或者其他特定值。

在此介紹一個新的套件：parameters，\(\tt{describe_distribution}\)可以列出連續變數的統計值以及區間估計，例如我們想知道所有大學收到的入學申請的平均份數加上\(95\%\)的區間估計，結果平均份數是3001，\(95\%\)的區間估計是(2786.31, 3195.65)。換句話說，假如某所學校收到的入學申請不到2786，那麼該學校就落在最後的5%的範圍。反過來假如某所學校收到的入學申請超過3195，那麼該學校就落在領先的5%的範圍。

 parameters::describe_distribution(
  ISLR::College$Apps,
  centrality = "mean",
  dispersion = FALSE, iqr = TRUE,
  range = FALSE, quartiles = FALSE,
  include_factors = FALSE, ci = 0.95, iterations = 100, threshold = 0.1 )
##    Mean |  IQR |             95% CI | Skewness | Kurtosis |   n | n_Missing
## ---------------------------------------------------------------------------
## 3001.64 | 2859 | [2699.84, 3240.38] |     3.72 |    26.77 | 777 |         0

• 我們模擬常態分佈的資料並且畫圖1.1呈現觀察值的分佈以及統計數字，再檢定平均數是否等於1：

library(ggstatsplot)
df <- tibble(X=rnorm(1000, 0.5, 1))
gghistostats(
  data = df,
  x = X, test.value = 1)

Figure 1.1: 長條圖與平均數統計

---

2 母體與樣本

研究某個現象時，我們假設這個現象來自於看不到的母體，要形容這個看不到的母體，我們會用這個母體的平均值以及離散程度這兩個參數。我們相信雖然我們永遠無法真正知道母體平均值以及離散程度，但是透過隨機抽樣的樣本，我們可以估計這兩個參數。

2.1 點估計步驟

點估計可按照以下步驟進行：

1. 抽取代表性的樣本
2. 選擇估計式(estimator)
3. 計算統計量的值
4. 以樣本統計量推論母體參數

2.2 估計式的選擇標準

如果母體參數是平均數，樣本的平均值是最好的估計方法嗎？樣本平均值是許多估計式(estimator)其中之一，但是我們為什麼要選擇平均值？估計式有以下的選擇標準：

1. 無偏估計：樣本統計量是否等於母體參數？是的話就是無偏估計。
2. 有效性：樣本統計量的平均平方誤差是否最小？也就是每次得到的估計值之間的離散程度最小。
3. 穩定(Consistency)：理論上當樣本數夠大時，估計值的分佈應該接近母體，如果不接近的話就不是穩定。
4. 常態(Asymptotic Normality): 當樣本數夠大時，統計樣本的分佈是否接近常態分佈？

• 偏誤是樣本估計\(\hat{\theta}\)與母體\(\theta\)之間的差距。

• 樣本估計的偏誤可表示為：\(\rm{Bias}(\hat{\theta})=E[\hat{\theta}-\theta]=E[\hat{\theta}]-\theta\)

• 無偏估計成立的條件是\(\rm{Bias}(\hat{\theta})=0\Longleftrightarrow \hat{\theta}=\theta\)

2.2.1 估計式的評估

• 根據以上四個原則可以評估什麼是好的統計樣本。不一定每個統計樣本都符合四種條件。

• 以下四種估計平均值的方式分別為:

1. 用其中一個樣本的平均數當做母體平均數，表示為：

\[\hat{\mu_{1}}=Y_{1}\]

2. 用一個常數當做母體平均數，表示為：

\[\hat{\mu_{2}}=c\]

3. 用樣本平均數的平均數做為母體平均數，表示為：

\[\hat{\mu_{3}}=\bar{Y_{n}}\equiv \frac{1}{n}(Y_{1}+Y_{2}+\cdots +Y_{N})=\mu\]

4. 分母+1的平均數做為母體平均數，表示為：

\[E[\hat{\mu_{4}}]\equiv \frac{1}{n+1}(Y_{1}+Y_{2}+\cdots +Y_{N})=\frac{n}{n+1}\mu \]

2.2.2 無偏估計

四種估計方式可分別計算平均數如下：

• \(E[\hat{\mu_{1}}]=E[Y_{1}]=\mu\)
• \(E[\hat{\mu_{2}}]=c\)
• \(E[\hat{\mu_{3}}]=\frac{1}{n}(E[Y_{1}]+E[Y_{2}]+\cdots +E[Y_{n}])=\frac{1}{n}(\mu+\mu+\cdots +\mu)=\mu\)
• \(E[\hat{\mu_{4}}]=E[\frac{1}{n+1}(Y_{1}+\cdots +Y_{n})]=\frac{1}{n+1}(E[Y_{1}]+E[Y_{2}]+\cdots +E[Y_{N}])=\frac{n}{n+1}\mu\)

\(\blacktriangleright\)由以上推導得知，只有\(\hat{\mu_{1}}\)與\(\hat{\mu_{3}}\)是無偏估計，也就是 \(\hat{\theta}=\theta\)。

2.2.3 有效性

在沒有偏誤的樣本估計中，選出一個\(\hat{\theta}\)有較小的變異的抽樣分布。

• \(\hat{\theta_{1}}\) 比 \(\hat{\theta_{2}}\)更穩定，如果 \(V(\hat{\theta_{1}})< V(\hat{\theta_{2}})\)

樣本統計的抽樣分佈的標準差稱為標準誤，寫成：

\[\sqrt{V[\hat{\theta}]}\]

回顧四種估計方式：

1. 用其中一個樣本的平均數當做母體平均數：\(\hat{\mu_{1}}=Y_{1}\)
2. 用一個常數當做母體平均數：\(\hat{\mu_{2}}=c\)
3. 用樣本平均數的平均數做為母體平均數：\(\hat{\mu_{3}}=\bar{Y_{n}}\equiv \frac{1}{n}(Y_{1}+Y_{2}+\cdots +Y_{N})\)
4. \(\hat{\mu_{4}}\equiv \frac{1}{n+1}(Y_{1}+Y_{2}+\cdots +Y_{N})\)

2.2.4 變異數

比較四種估計方式的變異數：

1. \(V[\mu_{1}]=V[Y_{1}]=\sigma^2\)（\(Y\)的變異數）
2. \(V[\hat{\mu_{2}}]=V[c]=0\)

\[\begin{align*} V[\hat{\mu_{3}}] & =V[\frac{1}{n}(Y_{1}+Y_{2}+\cdots +Y_{N})] \\ & =\frac{1}{n^2}(V[Y_{1}]+\cdots +V[Y_{n}]) \\ & =\frac{1}{n^2}(\sigma^2+\cdots +\sigma^2) \\ & =\frac{1}{n}\sigma^2 \end{align*}\]

4. \(V[\hat{\mu_{4}}]=\frac{n}{(n+1)^2}\sigma^2\)

\(\blacktriangleright\)比較\(\hat{\mu_{1}}\)以及\(\hat{\mu_{3}}\), \(\hat{\mu_{3}}\)的抽樣分佈變異數較小。

2.2.5 Mean Squared Error

在無偏估計中選擇一個穩定的樣本統計的標準是誤差平方的平均值(Mean Squared Error) ，或者稱為均方誤差。定義如下：

\[\begin{align*} MSE(\hat{\theta}) &=E[(\hat{\theta}-\theta)^2]\\ &=Bias(\hat{\theta})^2+V(\hat{\theta})\\ &=[E(\hat{\theta})-\theta]^2+V(\hat{\theta}) \end{align*}\]

• 均方誤差由偏誤（估計式的期望值與無法觀察的母數的差）以及估計式期望值的變異數所組成，其中一個大，另一個就會小。如果是無偏估計，偏誤等於0，那麼變異數等於均方誤差。

2.2.6 數理分析

1. \(MSE(\hat{\mu_{1}})=\sigma^2\)
   
2. \(MSE(\hat{\mu_{2}})=(c-\mu)^2\)

\[\begin{align*} MSE(\hat{\mu_{3}}) &=E[(\hat{\mu_{3}}-\mu)^2]\\ &=\rm{Bias}(\hat{\mu_{3}})^2+V(\mu_{3})\\ &=0+\frac{\sigma^2}{n} \end{align*}\]

4. \[\begin{align*} MSE(\hat{\mu_{4}})&=E[(\hat{\mu_{4}}-\mu)^2]\\ &=Bias(\hat{\mu_{4}})^2+V(\mu_{4})\\ &=\rm{Bias}(\frac{n}{n+1}\mu)^2+\frac{n}{(n+1)^2}\sigma^2\\ &=\frac{1}{(n+1)^2}\mu^2+\frac{n}{(n+1)^2}\sigma^2\\ &=\frac{\mu^2+n\sigma^2}{(n+1)^2} \end{align*}\]

• 以上數理分析的結論為：

\(\blacktriangleright\) \(\hat{\mu_{1}}\) 與 \(\hat{\mu_{3}}\) 都是無偏但是後者的誤差平方的平均值較小。換句話說，用樣本平均數的平均數做為母體平均數的估計式最好。

\(\blacktriangleright\) \(\hat{\mu_{3}}\)與\(\hat{\mu_{4}}\)很難說哪一個誤差平方的平均值較小，但是後者不是母體平均數的無偏估計。

2.2.7 樣本數趨近無限大的特性

當樣本數在一定的規模時，樣本統計的無偏以及穩定特性都成立。

當樣本數趨近無限大時，抽樣分佈會趨近於特定的區間，而樣本統計系列\(\theta_{1},\cdots \theta_{n}\) 則會趨近於特定的數，例如\(\mu\)。

2.2.8 穩定

穩定的樣本估計是\(\theta_{1},\cdots \theta_{n}\)會聚合在母體參數\(\theta\)，當\(n\rightarrow \infty\)，也就是：

• 理論上，可以分析\(E[\theta_{n}]\rightarrow \theta\)以及\(V[\theta_{n}]=0\)是否為真。

• 也可以用模擬的方式檢驗是否當樣本數增加，抽樣分佈是否趨近一直線。 例如圖2.1與圖2.2分別表示第三種與第四種樣本統計之模擬圖，樣本數=20以及樣本數=2000。

v1<-here::here('Fig','unbiasedness_vbig.jpg')
knitr::include_graphics(v1)

Figure 2.1: 樣本數20之模擬圖

v2<-here::here('Fig','unbiasedness_vsmall.jpg')
knitr::include_graphics(v2)

Figure 2.2: 樣本數2000之模擬圖

• 在樣本數趨近於無限大時，樣本統計的抽樣分佈會趨近於特定值。也就是變異數趨近於0。

• \(\hat{\mu_{3}}\)是無偏且穩定之樣本統計，也就是\(E[\hat{\mu_{3}}]=\mu\)。 而\(\hat{\mu_{4}}\)是偏差但穩定之樣本統計。

2.2.9 弱大數法則(Weak Law of Large Numbers)

弱大數法則指的是，如果有\(X_{1}, X_{2},\cdots X_{n}\)等一系列的隨機變數，彼此互相獨立，每一個變數的平均數等於母體參數\(\mu\)。對\(\epsilon>0\)，\(\bar{M_{n}}\)為平均數\(\frac{M_{1}+M_{2}+\cdots +M_{n}}{n}\)

• 當\(n\rightarrow \infty\)
• 以上的等式意義為當\(n\)趨近無限大，樣本平均數與母體平均數之間將沒有差異。
• 如果樣本數 \(n\) 夠大，那麼 \(M_{n}\) 的分佈將在 \(\mu\) 的附近，或者是 \([\mu-\epsilon,\mu+\epsilon]\)這個區間
• 當\(n\rightarrow \infty\)，樣本平均數會收斂於母體參數，也就是期望值。

2.2.10 強大數法則(Strong Law of Large Numbers, LLN)

強大數法則指的是，如果有\(X_{1}, X_{2},\cdots X_{n}\)等一系列的變數，彼此互相獨立，母體參數為\(\mu\)。對\(\epsilon>0\)，\(\bar{X_{n}}\)為平均數，則有以下結果：

• 或者是

\[\begin{align*} Pr({\mathrm {lim}}_{n\rightarrow \infty}|\bar{X_{n}}=\mu)=1 \end{align*}\]

• 強大數法則的意義為當\(n\rightarrow \infty\)，樣本平均數等於母體參數的機率會收斂於1。

接下來我們運用大數法則在母體平均數或者比例的估計。

2.3 母體平均數的估計

• 身高、體重、智商等變數的母體為常態分佈，平均數為\(\bar{X}=\frac{X_{1}+\ldots+X_{n}}{n}\)，機率密度函數為：

\[\begin{align*} f(X) & =\frac{1}{\sigma\sqrt{2\pi}}\cdot e^{\frac{-(x-\mu)^2}{2\sigma^2}} \end{align*}\]

• 常態分佈有以下性質：
• \(\bar{X}\sim N(\mu, \frac{\sigma^2}{n})\)
• \(\text{P}(\bar{X}-\mu \leq 1.96\sigma_{\bar{X}})=0.95\)

2.4 母體比例的估計

• 擲硬幣的結果、身分證最後一號是奇數或者偶數等等變數，屬於\(n\)次獨立重複伯努利試驗。事件\(X\)發生的次數表示為 \(n_{x}\)。\(\frac{n_{x}}{n}\)為事件\(X\)發生的機率，事件\(X\)在每次試驗中發生的母體機率為\(p\)。對任意正數\(\epsilon>0\)，以下等式成立：

\[\begin{align*} {\mathrm {lim}}_{n\rightarrow \infty}\text{Pr}(|\frac{n_{x}}{n}-p|< \epsilon)=1 \end{align*}\]

• 強大數法則的意義為當\(n\rightarrow \infty\)，樣本平均數等於母體參數的機率會收斂於1。請見圖2.3表示二元變數之抽樣分佈模擬結果。

v3<-here::here('Fig','unbiasbinomial41.jpg')
knitr::include_graphics(v3)

Figure 2.3: 二元變數之抽樣分佈模擬

• 丟硬幣100次，其中得到52次正面、48次反面。我們對於母體的平均數估計是\(\frac{52}{100}\cdot 1+\frac{48}{100}\cdot 0=0.52\)。

• 如果我們相信母體平均數是0.52，那麼擲這枚硬幣10次，其中會得到6次正面的機率是：\(\frac{10!}{6!4!}(0.52)^6\cdot (1-0.52)^4=\)

p=0.52; n=10; x=6
(factorial(n))/(factorial(x)*factorial(n-x))*p^x*(1-p)^4
## [1] 0.2204

• 假設我們抽樣600個同學，抽20次，而且我們相信女生佔了52%，那麼這50套樣本當中，有多少套的樣本，女生剛好52人？

set.seed(116)
m = 20; n= 600; p = .52;
phat = rbinom(m,n,p)/n
phat
##  [1] 0.5383 0.5200 0.5217 0.5250 0.4967 0.4833 0.5183 0.5233 0.4950 0.5517
## [11] 0.5467 0.5533 0.5400 0.5283 0.5300 0.5350 0.5350 0.5250 0.5067 0.5433

• 我們排序50次的結果，然後以點狀圖2.4表示：

dotchart(sort(phat), pch=16, col='#AE00AA', yaxt='n', xaxt='n',   
  ylab='',xlab=expression(paste("", hat(pi))), cex.lab=1.8)
abline(v = 0.52, lwd=1.5, col='red3', lty = 2)
#mtext('0.52',1)
axis(1, at=c( 0.49, 0.50, 0.51, 0.52, 0.53, 0.54,0.55),
     labels=c(0.49, 0.50, 0.51, 0.52, 0.53, 0.54,0.55))

Figure 2.4: 二項分布抽樣之點狀圖

• 圖 2.4 呈現出抽樣的結果，涵蓋小於0.45到大於0.6，但是只有幾個等於0.52。換句話說，抽樣得到的估計點，雖然理論上都是無偏估計，但是實際上與母體參數可能有一些不同。

2.5 中央極限定理(Central Limit Theorem, CLT)

統計學有兩大支柱，第一是大數法則，第二就是中央極限定理。大數法則告訴我們樣本平均數在無數次的實驗之後，會趨於根據期待值的線性法則(Linearity of Expectation)所得到的期待值與變異數，也就是\(\bar{X_{n}}\rightarrow \mu\quad \rm{as} \quad n\rightarrow \infty\)以及\(\rm{Var}(\bar{X_{n}})=\frac{\sigma^2}{n}\)。但是並沒有告訴我們期待值所形成的分佈。

根據中央極限定理，如果母體的分佈是常態，即使樣本數只有2個，樣本平均數的分佈會是常態分佈。如果母體的分佈不是常態，只要樣本數超過30個，樣本平均數的分佈會是常態分佈。如果母體的分佈不是常態，樣本數小於30個，樣本平均數的分佈會是\(t\)分佈。用圖2.5表示（圖的來源）：

include_graphics("https://miro.medium.com/max/700/1*qBenTmidG7bQdvsPdGjexA.png")

Figure 2.5: 中央極限定理示意圖

中央極限定理指的是：有\(X_{1}, X_{2},\cdots X_{n}\)等一系列的變數，彼此互相獨立，母體參數為\(\mu\), \(\sigma^2<\infty\)時，對任何隨機變數的母體分佈下列公式成立：

\[\begin{align*} \sqrt{n}(\bar{X_{n}}-\mu)\xrightarrow{d} N(0, \sigma^2) \end{align*}\]

• \(n\)變大，\(\sqrt{n}\)倍的樣本平均數將會聚合在常態分佈。
• 因為\(Var[ax]=a^2Var[x]\),也就是對一個變數乘上\(a\)倍,它的變異數便是乘上\(a^2\)倍。所以可得到(2.1)

\[\begin{align*} 1/\sigma \cdot \sqrt{n}(\bar{X_{n}}-\mu) \rightarrow 1/\sigma^2\cdot N(0, \sigma^2)=N(0,1)\tag{2.1} \end{align*}\]

• 接著，計算\(1/\sigma \cdot \sqrt{n}(\bar{X_{n}}-\mu)\)：

\[\begin{align} \frac{1}{\sigma}\cdot \sqrt{n}(\bar{X_{n}}-\mu) & = \frac{\sqrt{n}(\bar{X_{n}}-\mu)}{\sigma}\nonumber \\ & = \frac{1/\sqrt{n}\cdot \sqrt{n}(\bar{X_{n}}-\mu)}{1/\sqrt{n}\cdot \sigma}\nonumber \\ & =\frac{\bar{X_{n}}-\mu}{\frac{\sigma}{\sqrt{n}}}\tag{2.2} \end{align}\]

根據方程式(2.1)以及(2.2)：

\[\begin{align} Z_{n}\equiv \frac{\bar{X}_{n}-E[\bar{X}_{n}]}{\sqrt{V[\bar{X}_{n}]}} &=\frac{\bar{X}_{n}-\mu}{\sigma/\sqrt{n}}\xrightarrow{d}{\mathit N}(0,1)\tag{2.3} \end{align}\]

• 中央極限定理可以應用到樣本統計(\(\bar{X}_{n}\))。
• 方程式(2.3)顯示：當\(n\)到一定規模時，標準化常態分佈(Standardized Normal Distribution)成立

\[\begin{equation} \frac{|\bar{X}_{n}-\mu|}{\sigma/\sqrt{n}}\sim N(0,1) \end{equation}\]

• 意謂著

\[\begin{equation} \bar{X}_{n}\sim N(\mu,\sigma/\sqrt{n}) \end{equation}\]

• 當樣本規模夠大時，不論是來自於哪一種分佈，樣本統計會形成常態分佈的抽樣分佈。

• 常態分佈的抽樣分佈可以用標準常態分佈表示，\(\text{Z}\sim \text{N}(0, 1)\)。我們用圖 2.6 表示信賴水準0.95的區間：

Figure 2.6: 標準常態分佈信賴水準0.95的區間

• 圖 2.6 顯示每一個Z值對應的機率密度。我們可以畫\(\rm{Pr}(X\leq z)\)的累積機率圖如圖 2.7。

plot.new()
zx<-seq(-3, 3, length.out=100)
cat("P(Z=-3)=", pnorm(zx[1], 0, 1),"\n")
## P(Z=-3)= 0.00135
cat("P(Z=0)=", pnorm(zx[50], 0, 1),"\n")
## P(Z=0)= 0.4879
cat("P(Z=-3)=", pnorm(zx[100], 0, 1))
## P(Z=-3)= 0.9987
plot(zx, pnorm(zx, 0, 1), cex=0.3, col="#ea1021")
abline(h=pnorm(zx[50], 0, 1), v=0 )

Figure 2.7: Z值累積機率圖

• 二項分布指的是在\(n\)次的伯努力實驗中，如果用\(X\)隨機變數代表事件為真的次數，而且已知成功機率為\(p\)，當\(X=k\)的機率是：\({n}\choose{k}\) \(\times p^{k}(1-p)^{n-k}\)。我們用圖 2.8 表示當\(X\)從0到100而次數固定為100以及母體機率為0.5，對應的機率類似常態分佈，但是實際上這不是機率密度分佈，而是機率質量分佈，因為\(X\)並不是連續變數：

plot.new()
k <- c(0 : 100)
size=100; pb=0.5
y <- dbinom(k, size=size, prob=pb)
plot(k, y,  xlab = "count", ylab = "Pr(X=k)", type="n", main = "")
lines(k, y, col="#0011aa")

Figure 2.8: 二項分佈的PMF

• 圖 2.9 則是表示二項分佈\(\text{Pr}(X\leq k)\)的機率，也就是\(\sum_{x\leq k}\text{P}(X=k)\)的機率：

plot.new()
k <- c(0 : 100)
size=100; pb=0.5
y <- pbinom(k, size=size, prob=pb)
plot(k, y,  xlab = "count", ylab = "Pr(X=k)", type="n", main = "")
lines(k, y, col="#aa0022")

Figure 2.9: 二項分佈累積機率

2.6 信賴區間

雖然點估計是無偏估計，但是區間估計比點估計可靠。

區間估計需要考慮變數的平均數的變異程度，因為我們每次抽出的樣本所得到的平均數不會都相同，所以\(\bar{X}\)的變異數\(\sigma_{\bar{X}}\)的大小會影響區間估計。

區間的信賴水準(confidence level)以\(1-\alpha\)表示，\(\alpha=0.05\)時，代表信賴水準為0.95或者\(95\%\)。

因此，信賴區間是在一個既定信賴水準構成的區間，包含樣本統計量以及抽樣誤差。

步驟：

1. 選擇估計式、計算估計值
2. 取得樣本統計量的抽樣分配，例如常態分配、t分配。在常態分佈而且信賴水準為\(95\%\)時，母體參數\(\mu\)與樣本平均值之間的差距，應該不會大於\(1.96*\sigma_{\bar{X}}\)。
3. 得出母體參數的信賴區間
4. 得出母體參數的結論

我們用區間估計表示抽樣的結果，也就是點估計加上一定的區間。<ㄥli>

set.seed(116)
m = 50; n= 100; p = .52;
phat = rbinom(m,n,p)/n
SE = sqrt(phat*(1-phat)/n)
alpha = 0.05;zstar = qnorm(1-alpha/2)
cat("zstar=", zstar)
## zstar= 1.96
matplot(rbind(phat - zstar*SE, phat + zstar*SE),
  rbind(1:m,1:m), type="l", lty=1,
xlab=expression(paste(hat(p)-z,"*",SE,",  ", hat(p)+z,"*",SE)), ylab="")
abline(v=p) 

Figure 2.10: 區間估計圖

這個圖 2.10 是由以下幾個參數所組成：

• \(\hat{p}\) 表示每一套樣本當中，女生的比例
• SE表示樣本的標準誤=\(\sqrt{\frac{p(1-p)}{n}}\)
• \(\alpha\)(alpha)表示我們估計時容許的抽樣誤差，也稱為顯著水準(significance level)，顯著水準越小，代表誤判虛無假設成立的機會越小。一般設定為0.05或是\(5\%\)，分成左右兩邊各是0.025
• \(z^{*}\)表示一定的顯著水準下，對應的z值。顯著水準越大，z值越小。z值來自平均為0，變異數為1的標準常態分佈（iid \(N(0,1)\)）。
• \(z^{*}\cdot SE\)表示一定顯著水準下，涵蓋母體參數的上下限。顯著水準越小、範圍越大、區間越容易涵蓋母體參數。
• 從這個圖可以看出，大部分的線跨過0.52這條垂直線，代表大部分抽樣的結果，涵蓋了母體為0.52這個參數。
• 因為區間估計比點估計來得周全，所以我們在檢定平均數或者平均數差異時，將使用區間估計。
• 標準誤表示為\(\frac{\sigma}{\sqrt{n}}\)。連續變數的\(\sigma\)的無偏估計是\(s=\frac{\sum(X_{i}-\bar{X})^2}{n-1}\)。二元變數則是\(\sqrt{p(1-p)}\)。

整理以上的資訊可以得到：

• \(\hat{p}\)的抽樣分配寫成：\(\hat{p}\sim N(p, \frac{p(1-p)}{n})\)。當\(\hat{p}=0.5\)時，\(\frac{p(1-p)}{n}\)會最大，也就是抽樣誤差會最大，避免低估所以如果我們不知道母體比例\(p\)，可以先以\(\hat{p}=0.5\)代入。
• 利用\(Z\)分配求得\(\hat{p}\)與\(p\)的抽樣誤差為：\(\text{P}(|\hat{p}-p|\leq Z_{\alpha/2}\frac{p(1-p)}{n}) =1-\alpha\)
• 因此，母體比例\(p\)的信賴區間為：\([\hat{p}-Z_{\alpha/2}\frac{p(1-p)}{n}\leq p\leq \hat{p}+Z_{\alpha/2}\frac{p(1-p)}{n}]\)

2.6.1 母體平均數之區間估計（變異數已知）

在描述統計的課程中，我們介紹常態分佈可計算\(Z\)值，轉換成標準常態分佈之後，用\(Z\)值表示隨機變數的相對大小或者事件發生機率，例如薪水、身高等等。在計算\(Z\)值時，假設我們已經知道母體變異數\(\sigma\)，或者用樣本變異數代替。我們這裡分成已知母體變異數以及未知兩個部分來討論母體平均數的區間估計。

★統計一群受訪者的睡眠時間，得到數據如下：

• 睡眠時間的平均值為6.5小時

• 標準差: 0.5小時

• 樣本數: 100

• 請計算這一群受訪者的睡眠時間標準誤：

\[\sigma_{\bar{X}}=\frac{\sigma}{\sqrt{n}}=0.5/10=0.05\]

• 因為在常態分佈而且信賴水準\(95\%\)時，\(|\bar{X}-\mu|\leq 1.96*\sigma_{\bar{X}}\)，所以

\[|6.5-\mu|\leq 1.96*0.05=0.098\]

• 因為 \[\begin{equation*} |a-b|= \begin{cases} a-b & \text{if}\quad a-b\geq 0\\ b-a & \text{if}\quad a-b\leq 0 \end{cases} \end{equation*}\]

• 所以 \[\begin{equation*} |6.5-\mu|= \begin{cases} 6.5-\mu\leq 0.098 & 6.5-0.098=6.402\leq \mu \\ \mu-6.5\leq 0.098 & \mu \leq 6.5+0.098=6.598 \end{cases} \end{equation*}\]

• 因此，睡眠的平均時間應該落在\(6.402\leq \mu \leq 6.598\)的區間。

2.6.2 母體比例的區間估計（變異數未知）

★統計一群受訪者的政府滿意度，得到數據如下：

• 樣本數: 1013

• 滿意政府的人數為466

區間估計如下：

1. 計算樣本比例：\(\hat{p}=466/1013=0.46\)
2. 決定信賴水準：\(\alpha=0.05\)
3. \(z_{\alpha/2}=1.96\)

• 可用R計算如下：

alpha=0.05
zstar=-qnorm(alpha/2)
cat("z value=", zstar)
## z value= 1.96

4. 95%信賴區間：

\[\begin{align*} \hat{p}\pm 1.96\cdot \sqrt{\frac{0.46\cdot0.54}{1013}} & = 0.46\pm 1.96*0.015 \\ & = 0.46\pm 0.03 \end{align*}\]

5. 結論：母體比例區間估計為 \((0.46-0.03\leq p\leq 0.46+0.03)=(0.43\leq p\leq 0.49)\)

2.6.3 樣本數與信賴水準

• 因為 \(p=0.5\) 時，極大化標準誤 \(\sqrt{\frac{p\cdot (1-p)}{n}}=\sqrt{0.25}{n}=\frac{0.5}{\sqrt{n}}\) 。

• 如果信賴水準是0.95，\(z_{\alpha/2}\)經過查表或者計算，等於\(1.96\approx 2\)，抽樣誤差等於\(z_{\alpha/2}*\sqrt{\frac{p\cdot (1-p)}{n}}\approx 2\cdot \frac{0.5}{n}=\frac{1}{\sqrt{n}}\)。

• 因此，我們通常用\(\frac{1}{\sqrt{n}}\)決定抽樣誤差，反過來，我們也可以先決定抽樣誤差，再來決定樣本數。

2.7 t分佈區間估計

雖然\(s\)是標準差\(\sigma\)的無偏估計，但是\(s\)只是隨機亂數，並不是母體分佈。我們可以改用\(t\)分佈。\(t\)分布近似常態分佈，樣本數越大越接近常態分佈。

• 以下的圖2.11顯示\(t\)分佈的自由度為1, 5, 10, 30時的機率分佈，\(t\)分佈的自由度是\(n-1\)：

curve(dt(x, 30), from = -5, to = 5, col = "orange", 
      xlab = "t", ylab = "density", lwd = 2)
curve(dt(x, 10), from = -5, to = 5, col = "green2", add = TRUE, lwd = 2)
curve(dt(x, 5), from = -5, to = 5, col = "navyblue", add = TRUE, lwd = 2)
curve(dt(x, 1), from = -5, to = 5, col = "grey40", add = TRUE, lwd = 2)
legend("topleft", legend = paste0("DF = ", c(1, 5, 10, 30)),
       col = c("grey40", "navyblue", "green2", "orange"),
       lty = 1, lwd = 2)

Figure 2.11: 不同自由度的t分佈

• \(t\)分佈形狀取決於自由度，自由度代表已知統計值之後，觀察值可以變動的數目。\(t\)分佈的形狀代表它的離散程度，或者是樣本變異數，而變異數來自於樣本平均值(\(\bar{X}\))。當我們已知平均值，只要有一個樣本不變，其他\(n-1\)個觀察值改變，還是可以得到相同的平均值。所以自由度是\(n-1\)。

• 以下的式子表示母體參數的上下區間的機率為\(1-\alpha\)：

\[P(\bar{X}-t_{\alpha/2,n-1}S/\sqrt{n}<\mu<\bar{X}+t_{\alpha/2,n-1}S/\sqrt{n})=1-\alpha\]

• 三種信賴區間的\(t\)值表示如表 2.1：

Table 2.1: 三種信賴區間的t值

信賴區間	alpha	df.1	df.4	df.15	df.3000
90%	0.10	6.314	2.132	1.753	1.645
95%	0.05	12.710	2.776	2.131	1.960
99%	0.01	63.660	4.604	2.947	2.576

• R可以求得以上的機率，例如信賴區間為95%時，\(\alpha\)除以2，加上自由度：

qt(0.025, 3000)
## [1] -1.961

• \(t\)分佈的區間估計的步驟：

1. 計算\(T\)值
2. 計算\(p\)值
3. 或者比較\(t^{*}\)與\(T\)值

★ 假設10位受試者的體重資料如下，請問平均體重相當於173磅嗎？

• 先計算\(T\)：

weight <-c(175, 176, 173, 175, 174, 173, 173, 176, 173, 179)
n=length(weight)
x_bar=mean(weight)
s2=sum((weight-x_bar)^2)/(n-1)
c=173
T=(x_bar-c)/sqrt(s2/n)
print(T)
## [1] 2.762

• 已知\(T\)值，用pt(T, df)求p-value：

T=2.761805
pvalue=2*pt(T, df=9, lower.tail = FALSE)
pvalue
## [1] 0.02205

• 之所以pt(T,df)乘以2，是因為T對應的累積機率值是從最右邊一直累積過來，因此右邊的部分乘以2，代表左右兩邊各一半。

• 或者比較\(T\)與\(t^{*}\)，即可知道\(T\)是否落在拒斥域，並畫圖2.12：

T=2.761805
t.star=qt(0.975, 9)
t.star
## [1] 2.262

plot(seq(-4,4, length=200), seq(0,0.5, length=200), type='n',
     xlab='t', ylab='')
curve(dt(x, 10), from = -5, to = 5, col = "green2", add = TRUE, lwd = 2)
T=2.761805
t.star=qt(0.975, 9)
abline(v=T, lty=2, lwd=1.5)
abline(v=t.star, lty=3, lwd=1.5, col='#0044ab')
text(t.star-0.2, 0.15, expression(paste('t'^'*',"=2.262")))
text(T+0.2, 0.08, expression(paste('T=2.761')))

Figure 2.12: 比較T值

• 直接用t.test()進行\(t\)檢定：

t.test(weight, mu=173)
## 
##  One Sample t-test
## 
## data:  weight
## t = 2.8, df = 9, p-value = 0.02
## alternative hypothesis: true mean is not equal to 173
## 95 percent confidence interval:
##  173.3 176.1
## sample estimates:
## mean of x 
##     174.7

• 分析顯示p-value <0.05，也就是說我們可以拒斥\(\mu=173\)的假設。

• 我們用圖 2.13 表示信賴水準0.05的區間：

  

  Figure 2.13: 顯著水準0.05的區間

• 圖 2.13顯示，如果T=-1.98或者1.98，只有\(\frac{1-0.95}{2}=0.025\)的機率會發生。當T=2.76，右邊的區域剩下非常小的區域，換句話說，要觀察到T=2.76的機率非常小，如圖 2.14：

scale <- 0.1 
x <- seq(-4, 4, scale) 
df=100
y <- dt(x, df) 
plot(x, y, type = "l", main="t-Test, t = 2.76") 

linepos <- 2.76 
abline(v = linepos) 
cutpoint <- (max(x) - linepos) / scale 

xt <- x[(length(x)-cutpoint):length(x)] 
yt <- y[(length(y)-cutpoint):length(y)] 

# draw the polygon 

n <- length(xt) 
xt <- c(xt[1], xt, xt[n]) 
yt <- c(0,yt,0) 
polygon(xt, yt, col="red" ) 

Figure 2.14: T值等於2.76

• T值對應的機率在曲線底下的面積是1，因此T值對應的機率可以用累積機率的方式呈現如圖 2.15。當T=-3時，累積機率為0.0026或者0.26%，而到了T=0的累積機率則是0.5，T=3時，就是累積機率等於1。

tx<-seq(-3, 3, length.out=100)
cat("P(T=-3)=",pt(tx[1], df=30),"\n")
## P(T=-3)= 0.002695
cat("P(T=0)=",pt(tx[50], df=30),"\n")
## P(T=0)= 0.488
cat("P(T=-3)=", pt(tx[100], df=30))
## P(T=-3)= 0.9973
plot(tx, pt(tx, df=30), cex=0.3, col="sandybrown")

Figure 2.15: t分佈累積機率圖

3 假設檢定

假設檢定與區間估計的原理相同，但是假設檢定多了建立假設以及驗證假設的步驟，例如我們要檢定某個樣本估計值是否等於特定常數，步驟如下：

1. 假設\(\mu=c\)。

2. 決定顯著水準\(\alpha\)或者信賴區間\((1-\alpha)\%\)

3. 計算自由度\(n-1\)

4. 計算\(T=\frac{X-c}{s/\sqrt{n}}\)

5. 計算顯著水準\(\alpha\)與自由度相對應的\(t^{*}\)值。\(\alpha=5\%\)、自由度等於1,000時，約為1.962

qt(0.975, 1000)
## [1] 1.962

6. 當\(p(T)<\alpha\)，表示會發生c的機會非常小，可以拒斥\(\mu=c\)的假設。

7. 同樣的當\(T>t_{*}\)，表示會發生c的機會非常小，可以拒斥\(\mu=c\)的假設。

• 由以上步驟可知，假設檢定涉及到樣本估計值與欲檢驗值的差距，以及樣本估計的變異程度，還有信賴水準的大小。

3.1 雙尾檢定虛無假設與對立假設

• 虛無假設：假設母體參數或者其他參數為真，除非證明非真。寫成\(H_{\text{0}}\)。

• 例如\(H_{\text{0}}: \mu = 0\).

• 對立假設：相對於虛無假設，提出不同的假設，寫成\(H_{\text{1}}\)。例如\(H_{\text{1}}: \mu \neq 0\).

• 假設樣本來自於常態分佈，而且抽取時互相獨立，樣本估計值除以標準誤之後形成常態分佈，因此樣本估計值與虛無假設中的母體參數之間的差距，除以標準誤之後，我們可以檢驗\(T\)值是否大於樣本分佈的特定信賴水準的臨界點\(t^{*}\)。如果大於臨界點，表示要觀察到這樣的差距的機會很小，也就是我們可以說這樣的差距統計上很顯著。換句話說，我們可以下結論說樣本估計值並不會等於虛無假設中的參數，我們不接受虛無假設。

• 反過來說，如果小於臨界點，表示很容易觀察到這樣的差距，也就是兩者之間沒有差距的機會很大，所以結論是虛無假設成立。

• 圖 3.1 顯示接受域以及拒斥域。

Figure 3.1: 拒絕域與接受域

3.2 檢驗假設結果

兩個錯誤：根據樣本統計量來做推論時，可能會出現兩種錯誤，第一種是母體參數真的如同假設一樣，第二種則是母體參數跟假設不同：

• 型一錯誤(Type I error): 當\(H_{\text{0}}\)為真（與母體參數相同），但是拒絕\(H_{\rm{0}}\)。以\(\alpha\)表示型一錯誤發生的機率:

\[\alpha=\text{Pr}(\text{reject}\hspace{.2em}H_{\text{0}}|H_{\text{0}}\hspace{.2em} \text{is}\hspace{.2em} \text{true})\]

&#9733全體烏龍茶都是淨重150公克，假設每罐烏龍茶平均淨重150公克的情況下，\(H_{\rm{0}}: \mu \geq 150\)，抽取一組烏龍茶，秤重之後的平均淨重小於150公克，換句話說，\(\bar{X}\leq 150\)。如果根據這個結果，判斷\(H_{\rm{0}}\)不為真，也就是拒斥虛無假設，那麼我們就犯了型一錯誤。

• 型二錯誤(Type II error): 當\(H_{\text{0}}\)為偽（與母體參數不相同），但是沒有拒絕\(H_{\rm{0}}\)。以\(\beta\)表示型二錯誤發生的機率:

\[\beta=\text{Pr}(\text{not reject}\hspace{.2em}H_{\text{0}}|H_{\text{0}}\hspace{.2em} \text{is}\hspace{.2em} \text{false})\]

• 圖 3.2 顯示左側有\(\alpha=0.05\)的拒斥區域。當\(H_{\rm{0}}: \mu=\mu_{0}\)為真，我們卻拒絕這個假設，我們犯錯的機率為\(\alpha\)。而且當\(H_{\rm{0}}: \mu\geq \mu_{0}\)為真，我們卻拒絕這個假設，我們犯錯的機率也是\(\alpha\)。

Figure 3.2: 第一型錯誤

• \(\alpha\)越小代表我們容忍誤判\(H_{\rm{0}}\)不為真的程度越小，也就是唯有在很極端的情況下才會接受\(H_{\rm{1}}\)，不然大多數時候即使有可能犯錯卻寧可相信\(H_{\rm{0}}\)為真。

• \(\alpha\)也被稱為「顯著水準」。根據\(\alpha\)，我們決定臨界點，當樣本統計量介於負無限大以及臨界點或者介於臨界點與無限大之間，我們可以判斷\(H_{\rm{0}}\)不為真。

• 圖 3.3 顯示右側有\(\beta\)的拒斥區域。當\(H_{\rm{0}}: \mu=\mu_{0}\)不為真，我們卻接受這個假設，我們犯錯的機率為\(\beta\)。而且當\(H_{\text{0}}: \mu\leq \mu_{0}\)不為真，我們卻接受這個假設，我們犯錯的機率也是\(\beta\)。

Figure 3.3: 第二型錯誤

• 雙尾檢定的虛無假設是：

\[H_{\text{0}}: \mu = \mu_{0}\]

• 而對立假設則是：

\[H_{\text{1}}: \mu \neq \mu_{0}\]

3.2.1 雙尾檢定範例

• 目前台灣的女性受雇員工的全年總薪資中位數約為45.6萬。調查一家公司的薪資，樣本數為102人，其中年薪超過45.6萬的女性員工是54人，請問該公司女性員工年薪超過45.6萬的比例是否為\(p=0.5\)？

• 因為伯努利分佈的樣本變異數是\(p(1-p)\)，所以樣本標準誤的公式為：

\[SE=\sqrt{\hat{p}(1-\hat{p})/n}\]

• T值的公式為：

\[T=\frac{\hat{p}-c}{\sqrt{\hat{p}(1-\hat{p})/n}}\]

• 我們用R計算並且檢定：

n=102
p=54/102
c=0.5
T=(p-c)/sqrt(p*(1-p)/n)
T
## [1] 0.5951

• 輸入T值，計算對應的p值：

T=0.595; n=102
pvalue= pt(T, df=n-1, lower.tail = F)
pvalue
## [1] 0.2766

• 因為p值大於0.05，所以我們的樣本可能來自於\(p=0.5\)的母體。換句話說，我們無法拒斥\(H_{\rm{0}}:p=0.5\)的虛無假設。

v4<-here::here('Fig','t_alpha1.jpg')
knitr::include_graphics(v4)

Figure 3.4: 雙尾檢定圖

• 我們也可用 \(\tt{prop.test()}\) 計算區間以及檢定（\(\tt{t.test()}\) 需要原始資料）：

prop.test(54, 102, p=0.5, conf.level=0.95)
## 
##  1-sample proportions test with continuity correction
## 
## data:  54 out of 102, null probability 0.5
## X-squared = 0.25, df = 1, p-value = 0.6
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.4285 0.6281
## sample estimates:
##      p 
## 0.5294

• 因為p-value>0.05，因此我們可以接受\(c=0.5\)的假設。換句話說，這一套樣本來自於\(p=0.5\)的母體。

• 由於R的 $

• 輸入T值，改用標準常態分佈計算對應的p值，得到0.72：

T=0.595
pvalue= pnorm(T, 0, 1)
pvalue
## [1] 0.7241

• 因為當T=0.5時，\(p\)值約為0.69，所以當T接近0.6時，\(p\)值也隨之增加。
• 我們模擬常態分佈的資料，做雙尾檢定，如圖3.5，可以看出多數的值落在兩邊虛線的中間，

library(dplyr);library(ggplot2)
n=300; set.seed(02138)
upper=prop.test(160, n, p=0.5, conf.level=0.95)$conf.int[2]
lower=prop.test(160, n, p=0.5, conf.level=0.95)$conf.int[1]
tibble(x=rnorm(300, 0.54, 0.05)) %>%
  ggplot(aes(x=x)) +
  geom_histogram(fill='gray90', bins=30) +
  geom_density(col='#0011EE') +
  geom_vline(xintercept = c(lower, upper), lty=2, col="#EE00AA") +
  geom_vline(xintercept = 0.54, lty=1, col="#AA1100") +
  theme_bw()

Figure 3.5: 模擬資料之雙尾檢定

3.2.2 單尾檢定

有時候我們需要檢定樣本平均數是否大於或者小於某一個數。例如學生的平均身高是否高於160公分，或者BMI值是否低於24。這時候我們可以用單尾檢定。

單尾檢定需要把拒斥的區域集中在其中一邊，例如圖 3.6表示拒斥區域在左尾，這是我們假設母體參數大於假設值，除非我們發現用來推論母體的樣本統計值非常小，小到落在拒斥區域，那麼我們就不接受母體參數大於假設值，也就是母體參數其實小於假設值。如果顯著水準為0.05，那麼就是右邊0.05或是左邊0.05，而不是0.025。可以想像顯著水準0.05除以2所對應的\(t^{*}\) 要比0.05對應的\(t^{*}\)來的大。

qnorm(0.025, 0, 1, lower.tail = T)
## [1] 1.96
qnorm(0.050, 0, 1, lower.tail = T)
## [1] 1.645

• 左尾檢定的虛無假設是：

\[H_{\text{0}}: \mu\geq \mu_{0}\]

• 而對立假設則是：

\[\text{H}_{\text{1}}: \mu < \mu_{0}\]

• 也可以詮釋成當我們用平均值\(\mu_{0}\)來隨機抽樣，要抽到\(<\mu_{0}\)（左尾）或是\(>\mu_{0}\)（右尾）的樣本的機會是\(p\)值。如果\(p\)小於5%，那麼我們可以說機會非常小，無法接受虛無假設成立。所以結論是我們可以接受對立假設，也就是\(\mu>\mu_{0}\)（左尾）或是\(\mu < \mu_{0}\)（右尾）。

• 例如過去研究宣稱某種病毒在人體表面可以生存超過5天，不管是6天、10天都比5天來得長。雖然我們懷疑這個研究發現過度誇大病毒的存活能力，但是我們除非發現的樣本統計值，換算成Z值或者t值之後小於臨界值，那麼我們才會推翻虛無假設。但是如果Z值或者t值大於臨界值，不論有多大，都符合虛無假設。

• 通常我們設定的標準是5%，要小於這個標準我們才拒斥虛無假設，也就是我們已經有點相信虛無假設，除非我們發現很有力的證據才能推翻虛無假設。圖3.6顯示，如果\(Z\)值小於-1.96，那麼我們可以拒斥\(H_{\text{0}}: \mu\geq \mu_{0}\)這個假設。

plot.new()
x <- seq(-4, 4, length = 1000)
y <- dnorm(x, 0, 1)
plot(x, y, type="n", xlab = "", ylab = "", main = "", axes = FALSE)
axis(1)
lines(x, y)

alpha=0.05
lower_bound <- 0-qnorm(1-alpha, 0, 1)
upper_bound <- 0+qnorm(1-alpha, 0, 1)
bounds_filter <- x >= -4 & x <= lower_bound
x_within_bounds <- x[bounds_filter]
y_within_bounds <- y[bounds_filter]
 
x_polygon <- c(-4, x_within_bounds, lower_bound)
y_polygon <- c(0, y_within_bounds, 0)
 
polygon(x_polygon, y_polygon, col = "red")
text(-2.6, 0.1, expression(paste(bold(alpha))))
segments(-2.6, 0.09, -2.2, 0.02, lty=2)
segments(0, 0, 0, 0.4, lty=2)

Figure 3.6: 左尾檢定圖

• 反之，如果我們想要檢定平均數小於\(\mu_{0}\)，那麼可以採用右尾檢定，也就是如果我們觀察到的樣本估計值非常大，大到母體平均數很難超過它，那麼我們可以拒斥平均值小於\(\mu_{0}\)的虛無假設，也就是接受大於\(\mu_{0}\)的對立假設。

• 右尾檢定的虛無假設是：

\[H_{\text{0}}: \mu\leq \mu_{0}\]

• 而對立假設則是：

\[H_{\text{1}}: \mu > \mu_{0}\]

• 圖3.7表示右側的區域為0.05。

v5<-here::here('Fig','t_alpha2.jpg')
knitr::include_graphics(v5)

Figure 3.7: 右尾檢定圖

• 上面的圖顯示當右尾的拒斥區域等於 5%，\(t_{\alpha}\)=1.69。

alpha=0.05
qt(1-alpha, df=29)
## [1] 1.699

3.2.3 單尾檢定範例1

• 在 UsingR::babies 這筆資料中，有父親的身高資料。請問父親的平均身高至少70.5英吋（約179公分）嗎？還是隨機誤差？請先扣掉99這個遺漏值。

• 因為我們感興趣的是父親身高是否高於70.5英吋，至於有多高並不是重點，所以我們要進行左尾的單尾檢定，對立假設是\(\mu <70.5\)，虛無假設則是：

\[H_{\text{0}}: \mu \geq 70.5\]

• 對立假設是

\[H_{\text{1}}:\mu < 70.5 \]

• 用R進行左尾的單尾檢定假設時，設定\(\mu\)等於\(\mu_{0}\)，而且\(\tt{alternative='less'}\)。

x <- UsingR::babies$dht
x <- x[x<99]
t.test(x, mu=70.5, alternative="less")
## 
##  One Sample t-test
## 
## data:  x
## t = -2.8, df = 743, p-value = 0.003
## alternative hypothesis: true mean is less than 70.5
## 95 percent confidence interval:
##   -Inf 70.38
## sample estimates:
## mean of x 
##      70.2

• 結果顯示，平均值為70.2，95%信賴區間的臨界值是70.37，小於70.5英吋，\(p\)值小於0.05，也就是說抽到的樣本其平均值會大於70.5的機會非常小。換句話說，我們可以拒絕虛無假設，並且結論父親的平均身高不到70.5英吋。

• 我們也可以手動計算

1. 以\(\bar{X}\)為檢定統計量，上面已知\(\bar{X}=70.2043\)。
2. 計算\(Z\)值

\[Z=\frac{\bar{X}-\mu_{0}}{S/\sqrt{n}}\]

x <- UsingR::babies$dht
x <- x[x<99]
n=length(x); cat("n=",n, "\n")
## n= 744
cat('SD=', sd(x), '\n')
## SD= 2.891
mean.x=mean(x); cat('Mean=',mean.x, '\n')
## Mean= 70.2
cat('Z=', (mean.x-70.5)/(sd(x)/sqrt(n)))
## Z= -2.79

• 計算出來\(Z\)值為-2.789，也就是落在左端與\(Z=0\)的中間。所以接下來用\(\tt{pnorm(Z, lower.tail=TRUE)}\)計算\(p\)值。

Z = -2.789
pnorm(Z, lower.tail = TRUE)
## [1] 0.002644

• 結果顯示，\(p < 0.05\)，表示可以拒斥虛無假設，也就是母體參數不會至少大於70.5，反之，應該是小於70.5。我們用機率密度圖3.8表示資料的分佈、檢驗值（70.5）以及左尾的臨界點(70.3)。

dht <- UsingR::babies$dht
dht <- dht[dht<99]
mu=70.5
Test<-t.test(dht, mu = mu, alternative="less")
upper=Test$conf.int[2]
lower=Test$conf.int[1]
upper; lower
## [1] 70.38
## [1] -Inf
tibble(dht) %>%
  ggplot(aes(x=dht)) +
  geom_density(col='#0011EE') +
  geom_vline(xintercept = upper, lty=2, col="#0b1abb") +
 # geom_vline(xintercept = mean(dht), lty=2, col='#AA1100') +
  geom_vline(xintercept = mu, col='#1100BB') +
  theme_bw()

Figure 3.8: 父親身高左尾檢定圖

• 我們也可以用\(t\)分佈檢定，也就是用\(\tt{pt()}\)，得到的結論相同，但是\(\tt{pt()}\)需要輸入自由度。

pt(-2.789, df=743)
## [1] 0.002711

☛請檢定\(\mu\geq70\)。

3.2.4 單尾檢定範例2

• 某家飲料工廠想要知道裝瓶的飲料是否符合包裝所標示的500ml，抽出30瓶測量其體積如下（參考Sarah Stowell）：

volume<-c(494.11,459.49,491.38,512.01,494.48,511.63,
  513.64,495.03,483.88,503.59,512.85,508.08,485.68,
  495.03,475.32,519.41,510.13,494.32,499.08,489.27,
  516.64,485.03,479.18,509.59,512.85,518.08,465.68,
  491.49,502.38,519.01)
length(volume)
## [1] 30

• 因為我們感興趣的是樣本平均值是否小於500ml，也就是飲料只能多不能少，所以進行右尾的單尾檢定。虛無假設是樣本平均值應該等於或小於500ml。

\[H_{0}: \mu\leq 500\]

\[H_{1}: \mu > 500\]

• 以R進行\(t\)檢定

t.test(volume, mu=500, alternative="greater")
## 
##  One Sample t-test
## 
## data:  volume
## t = -0.59, df = 29, p-value = 0.7
## alternative hypothesis: true mean is greater than 500
## 95 percent confidence interval:
##  493.4   Inf
## sample estimates:
## mean of x 
##     498.3
t.test(volume, mu=500)
## 
##  One Sample t-test
## 
## data:  volume
## t = -0.59, df = 29, p-value = 0.6
## alternative hypothesis: true mean is not equal to 500
## 95 percent confidence interval:
##  492.4 504.2
## sample estimates:
## mean of x 
##     498.3

• 結果顯示，\(t\)值為-0.594，\(p>0.05\)，也就是無法拒斥虛無假設。95%的信賴水準構成的區間顯示，母體平均值最小為493.36，樣本平均值為498.278，因此母體平均值可能小於500，但是應該會大於493。換句話說，我們無法拒斥樣本平均值應該等於或小於500ml的虛無假設。

• 用\(\tt{pt(T, df, lower.tail=FALSE)}\)計算\(p\)值，得到0.721，也就是右尾檢定的\(p\)值。

T=-0.59499
pt(T, 29, lower.tail =FALSE)
## [1] 0.7218

• 圖 3.9 呈現雙尾檢定的結果，其中\(p=0.556\)，這是因為進行雙尾檢定時，\(p\)值=\(2\times \tt{pt(T, df, lower.tail=FALSE)}\)。而95%的區間估計為(492.19, 504.37)，包含了樣本平均值498.28以及500，顯示母體參數可能大於或者小於500。

# for reproducibility
set.seed(123)

# plot
tibble(volume) %>%
ggstatsplot::gghistostats(
  data = ., 
  x = volume, 
  title = "", 
 # caption = substitute(paste(italic("Source:"), "")),
  type = "robust", # one sample t-test,
  test.value = 500, # default value is 0
  tr = 0,
  #conf.level = 0.95,
  bf.message = FALSE, # display Bayes method results
 # test.value.line = TRUE, 
  centrality.parameter = "mean", 
  centrality.line.args = list(color = "darkred"), 
  binwidth = 0.8, # binwidth value (experiment)
  messages = FALSE, # turn off the messages
  ggtheme = ggthemes::theme_tufte(), 
  ggstatsplot.layer = T,
  bar.fill = "#D55E00", # change fill color
)

Figure 3.9: 結合t檢定之長條圖

3.3 小樣本的平均數檢定

• 如果樣本小於30時，必須使用\(t\)檢定。步驟與上述的雙尾檢定以及單尾檢定一樣。

3.4 母體比例的檢定

• 母體參數可能是比例\(p\)，例如性別的比例、產品的不良率等等。跟區間估計一樣，檢定統計量\(Z\)或者\(T\)需要用到樣本的比例\(\hat{p}\)，也就是：

\[Z=\frac{\hat{p}-p}{\text{SE}(p)}\]

• 其中

\[\text{SE}(p)=\sqrt{\frac{p\times(1-p)}{n}}\]

1. 選定樣本分配：

\[\hat{p}\sim N(p,\frac{p(1-p)}{n})\]

2. 計算標準誤
3. 計算檢定值
4. 如果檢定值>0.05，則接受虛無假設。反之，則拒斥虛無假設。

3.4.1 範例

• 隨機抽樣1,000位大學生，其中615人表示希望繼續就讀研究所，請問所有學生就讀研究所的意願是否超過6成？

• 虛無假設：

\[H_{0}: P\leq 0.6\]

• 對立假設：

\[H_{1}: P>0.6\]

• 因為對立假設是\(P>0.6\)，所以要進行右尾檢定。我們先計算\(\hat{p}\)，然後計算標準誤，接著計算\(T\)值，並且計算\(p\)值。

phat=615/1000
p=0.6
n=1000
se=sqrt(p*(1-p)/n)
T=(phat-0.6)/se
cat("T=",T)
## T= 0.9682
cat("P=",pt(T, n-1, lower.tail = F))
## P= 0.1666

4 指令整理

Scott Stevens在他的書Introduction to Statistics: Think & Do整理了假設檢定的指令，請見Hypothesis Testing。

• \(\tt{qnorm(conf, mean, sd)}\): 已知平均數、變異數，求一定信賴水準下的Z值(\(\rm{z}_{\alpha/2}\))

alpha=0.01; cat("conf=0.99, mean=0, sd=1, Z=", qnorm(1-(alpha/2), 0, 1),"\n")
## conf=0.99, mean=0, sd=1, Z= 2.576
alpha=0.05; cat("conf=0.95, mean=0, sd=1, Z=", qnorm(1-alpha/2),"\n")
## conf=0.95, mean=0, sd=1, Z= 1.96
alpha=0.1; cat("conf=0.90, mean=0, sd=1, Z=", qnorm(1-alpha/2),"\n")
## conf=0.90, mean=0, sd=1, Z= 1.645
alpha=0.05; cat("conf=0.95, mean=1, sd=1, Z=", qnorm(1-alpha/2, 1, 1),"\n")
## conf=0.95, mean=1, sd=1, Z= 2.96

• 請參考圖 4.1 顯示信賴區間為0.95(\(\alpha=0.05\))的常態分佈，以及兩側各為0.025的拒斥區域。

Figure 4.1: 信賴區間0.95常態分佈圖

• \(\tt{pnorm(z, mean, sd, lower.tail=TRUE)}\): 已知平均數、變異數，求小於\(z^{*}\)的累積機率值(\(\text{p}(Z\leq z^{*})\))，參考圖 4.2。

cat("z=1.99, mean=0, sd=1, p=", pnorm(1.99, 0, 1, lower.tail = TRUE),"\n")
## z=1.99, mean=0, sd=1, p= 0.9767
cat("z=1.96, mean=0, sd=1, p=", pnorm(1.96, 0, 1, lower.tail = TRUE),"\n")
## z=1.96, mean=0, sd=1, p= 0.975
cat("z=1.68, mean=0, sd=1, p=", pnorm(1.68, 0, 1, lower.tail = TRUE),"\n")
## z=1.68, mean=0, sd=1, p= 0.9535

Figure 4.2: 左尾累積機率圖

• \(\tt{pnorm(z, mean, sd, lower.tail=FALSE)}\): 已知\(z^{*}\)值，平均數、變異數，求右尾機率值(\(\text{p}(Z\geq z^{*})\))，參考圖 4.3。

cat("z=1.99, mean=0, sd=1, p=", pnorm(1.99, 0, 1, lower.tail = FALSE),"\n")
## z=1.99, mean=0, sd=1, p= 0.0233
cat("z=1.96, mean=0, sd=1, p=", pnorm(1.96, 0, 1, lower.tail = FALSE),"\n")
## z=1.96, mean=0, sd=1, p= 0.025
cat("z=1.68, mean=0, sd=1, p=", pnorm(1.68, 0, 1, lower.tail = FALSE),"\n")
## z=1.68, mean=0, sd=1, p= 0.04648

Figure 4.3: 右尾累積機率圖

• \(\rm{qt(1-alpha/2, df)}\): 已知自由度，求一定信賴水準下的T值(\(\text{t}_{\alpha/2}\))

n=10; alpha=0.01; cat("conf=0.99, n=10, T=", qt(1-alpha/2, n-1),"\n")
## conf=0.99, n=10, T= 3.25
n=10; alpha=0.05; cat("conf=0.95, n=10, T=", qt(1-alpha/2, n-1),"\n")
## conf=0.95, n=10, T= 2.262
n=10; alpha=0.1; cat("conf=0.90,  n=10, T=", qt(1-alpha/2, n-1),"\n")
## conf=0.90,  n=10, T= 1.833
n=30; alpha=0.05; cat("conf=0.95, n=30, T=", qt(1-alpha/2, n-1),"\n")
## conf=0.95, n=30, T= 2.045

• \(\tt{pt(t, df)}\): 已知自由度，求累積機率值(\(\rm{p}(T\leq t^{*})\))

n=10; alpha=0.01;cat("t=1.99, n=10, p=", pt(1.99, n-1),"\n")
## t=1.99, n=10, p= 0.9611
n=10; alpha=0.05;cat("t=1.96, n=10, p=", pt(1.96, n-1),"\n")
## t=1.96, n=10, p= 0.9592
n=10; alpha=0.1;cat("t=1.68, n=10, p=", pt(1.68, n-1),"\n")
## t=1.68, n=10, p= 0.9364
n=30; alpha=0.05;cat("t=1.68, n=30, p=", pt(1.68, n-1),"\n")
## t=1.68, n=30, p= 0.9481

• \(\tt{pt(t, df)}\): 已知自由度，求右尾機率值(\(\text{p}(T\geq t^{*})\))

n=10; alpha=0.01;cat("t=1.99, n=10, p=", pt(1.99, n-1, lower.tail = FALSE),"\n")
## t=1.99, n=10, p= 0.0389
n=10; alpha=0.05;cat("t=1.96, n=10, p=", pt(1.96, n-1, lower.tail = FALSE),"\n")
## t=1.96, n=10, p= 0.04082
n=10; alpha=0.1;cat("t=1.68, n=10, p=", pt(1.68, n-1, lower.tail = FALSE),"\n")
## t=1.68, n=10, p= 0.06363
n=30; alpha=0.05;cat("t=1.68, n=30, p=", pt(1.68, n-1, lower.tail = FALSE),"\n")
## t=1.68, n=30, p= 0.05185

• \(\tt{t.test()}\): 針對原始資料，計算一定信賴水準下的信賴區間

set.seed(02138)
X<-rnorm(100, 0, 1)
t.test(X, conf.level=0.95)
## 
##  One Sample t-test
## 
## data:  X
## t = -0.81, df = 99, p-value = 0.4
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.2795  0.1168
## sample estimates:
## mean of x 
##  -0.08133
#manual computing
n=length(X)
mean.x<-mean(X)
SD <- sqrt(var(X)/n)
alpha=0.05
zstar<-qnorm(1-alpha/2)

cat("conf=0.95, interval estimate=[", mean.x-zstar*SD, ",", 
         mean.x+zstar*SD,"]")
## conf=0.95, interval estimate=[ -0.277 , 0.1144 ]

• \(\tt{prop.test()}\): 針對有比例的原始資料，計算一定信賴水準下的信賴區間(\(\hat{p}\pm t_{\alpha/2}\cdot \sqrt{\frac{p\times(1-p)}{n}}\))

set.seed(02138)
X<-rbinom(100, 1, 0.5)#randomly drawing samples from binomial distribution
n=length(X)
s<-length(X[X==1])
prop.test(s, n, conf.level = 0.95)
## 
##  1-sample proportions test with continuity correction
## 
## data:  s out of n, null probability 0.5
## X-squared = 0.49, df = 1, p-value = 0.5
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
##  0.3609 0.5622
## sample estimates:
##    p 
## 0.46
#manual computing
phat<-mean(X)
SE <- sqrt(phat*(1-phat)/n)
alpha=0.05
tstar<-qt(1-alpha/2, df=n-1)

cat("conf=0.95, interval estimate=[", phat-tstar*SE, ",", 
         phat+tstar*SE,"]")
## conf=0.95, interval estimate=[ 0.3611 , 0.5589 ]

• \(\tt{binom.test()}\): 針對有成功數的二元資料，計算一定信賴水準下的信賴區間(\(\hat{p}\pm t_{\alpha/2}\cdot \sqrt{\frac{p\times(1-p)}{n}}\))，結果與\(\tt{prop.test()}\)相近。

set.seed(02138)
X<-rbinom(100, 1, 0.5)#randomly drawing samples from binomial distribution
n=length(X)
s<-length(X[X==1])
binom.test(s, n, conf.level = 0.95)
## 
##  Exact binomial test
## 
## data:  s and n
## number of successes = 46, number of trials = 100, p-value = 0.5
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.3598 0.5626
## sample estimates:
## probability of success 
##                   0.46

5 小結

區間估計與假設檢定的原理相同，只是假設檢定需要虛無假設。

假設檢定又可分雙尾、左尾、右尾等三種。

雙尾檢定要用\(z_{\alpha/2}\)或者\(t_{\alpha/2}\)做為臨界值，左尾則是\(z_{\alpha}\)，右尾則為\(z_{1-\alpha}\)，其中\(\alpha=0.05\)。

6 作業

6.1 區間估計

1. 有一份調查的樣本數是1,000，性別的比例\(\hat{p}=0.45\)。請列出信賴水準為50%, 60%, 70%, 80%, 95%的區間估計。

2. 環保署監測某個縣市的二氧化碳每小時濃度，平均值為0.4，變異數為0.16。如果二氧化碳每小時濃度成常態分佈，請問前10%的每小時二氧化碳濃度至少是多少？

3. 隨機抽樣調查 350位同學，發現140位受訪者在最近一年曾經出國。請問信賴水準90%的區間估計是多少?

4. 有兩家民調公司，分別接受同一家廠商委託調查對於最新銷售的電腦的好感度。第一家公司收集\(n_{1}\)的隨機樣本，樣本變異數為\(\hat{\sigma_{1}}\)。第二家公司則認為只需要一半的隨機樣本即可，而且樣本變異數是另一家公司的四分之一，也就是\(\hat{\sigma_{2}}=0.25\hat{\sigma_{1}}\)。請回答以下子題：

1. 假設第一家公司的樣本數為60，調查得到的平均為\(\hat{\mu}_{1}=34\)，第二家公司調查得到的支持度為\(\hat{\mu}_{2}=39\)。請問這兩家公司對於產品好感度的最佳估計是什麼？

2. 請問這兩家調查的支持度信賴區間分別為何？

5. 請畫圖檢視\(\textbf{UsingR::brightness}\) 的分佈，並且以 \(\tt{t.test}\) 計算95%信賴水準的區間估計。

7. 請回答現代統計學二版的第9.10題。

8. alr4::walleye有walleye（玻璃梭鱸）的身長資料(length)，請問95%的區間估計為何？

9. 請問\(\text{P}(0<Z<1)=?\)

10. 隨機抽查1154位民眾，有48%的受訪者表示應該採用電子投票，有40%的受訪者認為應該維持紙本投票，另外有12%沒有意見。針對電子投票的贊成比例，請問信賴水準99%的區間估計為何？

6.2 假設檢定

1. 請檢定UsingR中的rat這筆資料的平均數是否大於110。

2. 請問UsingR中的<span style:“color=blue”>iq這筆資料中，iq的平均數大於100嗎？

3. 請回答現代統計學二版的第10.8題。

4. 請回答現代統計學二版的第10.12題。

5. 請回答現代統計學第10.16題第1小題。

6. 請回答現代統計學二版的第10.18題。

7. 請寫出右尾、左尾、雙尾等三種檢定的\(Z\)與\(Z_{\alpha}\)之間的關係。

8. 某社區醫院統計其接生的新生兒平均重量為2.8千克。若隨機選擇最近一個月在該醫院出生的10個嬰兒作為樣本，記錄他們的體重分別為2.8，4.0，3.3，4.1，2.9，2.8，2.9，2.7，3.3，3.2千克。請問如果顯著水準為5%，可否結論最近一個月出生的嬰兒平均體重等於2千8百公克？

9.承上題，如果有位研究人員看到前面幾個嬰兒的體重，便假設這幾個新生兒的平均體重一定超過3千6百公克，請問她的假設正確嗎？

10. 某電商平台規定，有問題的貨品比例不能超出7%，而倉儲抽出360個商品調查，顯示有5%的貨品有問題。倉儲懷疑這整批貨的NG比例超過7%。試問這批貨需要被退貨嗎？

7 更新內容日期

## 最後更新日期 05/06/2021