1 母體與樣本
- 1.1 抽樣架構(sampling frame)
- 1.2 無反應
2 描述統計
3 離散程度
4 標準常態分佈的Z值
- 4.1 Z值或分數
  - 4.1.1 範例一
  - 4.1.2 範例二
5 作業
6 更新講義時間

1 母體與樣本

國家用統計來了解全國的生產力以及土地、森林等資源，所以每一個單位都要納入範圍中。而社會科學則假設人類行為跟態度可以被測量，並且假設這些觀察值來自於某種分布，例如常態分佈、均等分佈、博桑分佈等等，然後用樣本統計來推論全體，並且找出哪些變數之間有相關。

例如我們想研究台灣有多少民眾了解中國，可以從電話簿、地址等等抽樣出一群民眾，然後調查他們的態度，記錄他們的反應，再來推論全體的民眾。態度可以視為連續變數，來自於常態分佈。

如果我們也想研究哪些國家發生民主倒退，假設倒退是0或者1的二元變數，也就是來自二元分佈，我們可以計算民主倒退的比例以及離散程度，然後找出跟這個現象相關的變數。
又例如我們想研究台北市的便利商店在午夜12點到凌晨5點之間，有多少人進入便利商店，我們可以抽樣100家超商，然後統計人數，推論到全部的超商，並且用超商的各種位置特徵，解釋這些超商樣本的人數，然後用來預測其他超商在深夜的來客數。人數屬於計數(count)資料，需要用博桑分佈來計算其平均數與離散程度。
有時候我們可以訪問母體裡的所有個人，這時候我們就不需要樣本。例如我們研究的是一本書裡面所有的文字背後的標題，只要我們有電子檔案，就可以分析所有的文字，而標題應該來自於Latent Dirichlet Allocation (LDA)或者多元分佈(multivariate distribution)。每一個標題都是一個選項，加總這些選項發生的機率等於1。

通常我們只會從母體抽出一次的樣本，樣本可以是一個也可以是幾百個或者上千個，樣本的特徵理論上不可能完全等於母體，也就是兩者之間會有誤差。但是我們如果抽無數次樣，每一套樣本都進行統計，例如計算平均值，那麼我們就會得到無數個平均值。理論上，樣本平均值的平均值將會等於母體的平均值。可以用以下式子表示：

\[ E[Y]=E[E(\hat{Y})] \]

1.1 抽樣架構(sampling frame)

抽樣需要一套規則，例如隨機抽樣，把所有人編號之後，隨機抽出若干號碼，然後進行調查。隨機抽樣的最大優點是理論上我們可以根據這個規則重複抽樣，每套樣本的平均數的平均數應該等於母體的平均數。

另外一個方法是設定一個樣本架構，按照規則抽樣。例如為了瞭解本校學生對於未來生涯的選擇，對於本校的大學生進行調查，可以先按照本校的學院進行分層，例如商學院、社科院的學生人數較多，可能要抽樣比較多的學生，但是外語學院、教育學院的學生人數比較少，抽樣的人數就比較少。在這樣的抽樣架構底下，再根據系級分層，或者隨機抽樣，大的學院分配到比較多的樣本，但是每個學生的中選機率比較低，反過來，小的學院的樣本數少，但是中選機率高，讓全校每一個學生的中選機率一樣，就可以兼顧隨機抽樣以及母體結構的需求。如果完全隨機抽樣，樣本可能散佈在每個學院，結果大的學院不見有比較多的樣本，雖然可以代表全校學生，但是沒有反應學院的差異，而不同學院的學生可能對於生涯發展有不同的想法。

第三種方法是集群(cluster)抽樣。例如市政府想了解文山區民眾對於文山區運動中心的滿意度，而且問卷的其中一個議題是使用文山運動中心的程度。為了這個目的，也許可以隨機抽台北市的電話簿，然後打電話詢問是不是住在文山區，是的話再詢問滿不滿意文山區運動中心。這樣做的一個可能缺點是要花很多時間篩選合格的受訪者。如果我們把文山區劃分為若干比較小的地理區，然後每個地理區裡面抽出若干家戶進行訪問，這樣做的好處是可以確實了解不同居住地點的民眾去運動中心的頻率，而且每一區都有一些樣本數。缺點是成本比高，而且要考慮不同地理區之間的差異。

還有一種抽樣方法是配額抽樣(quota sampling)。假設全台灣合格選民1900萬之中，男女各一半，20到29歲佔了15%，也就是說大概有142萬的20到29歲的男性，佔了全體合格選民的7.5%。假設我們想訪問1,000位的合格選民，那麼根據配額抽樣，我們應該訪問1000*7.5%=75位20到29歲的男性。理論上，我們可以在訪問結束時就得到符合母體人口比例的樣本。

配額抽樣的缺點是在缺乏適當的監督下，可能有因為訪員造成的選樣偏誤。如果是面訪，訪員自己根據配額去接觸受訪者，有可能遇到屬於這個配額的受訪者，但是因為該受訪者不太合作，所以就更換下一個屬於該配額的受訪者，造成非機率抽樣。因此，我們沒辦法反覆地根據這套方式反覆抽樣。不過，很多網路調查進行配額抽樣，如果已經確定某個性別、年齡層、教育程度以及其他人口特徵所佔的比例，而且有這些受訪者的資料庫夠大，的確可以透過網路邀請資料庫中的受訪者來填答，一直到配額額滿。

1.2 無反應

訪問過程中有可能遇到不合作的受訪者，拒絕接電話或者拒絕開門讓訪員訪問，我們稱為unit nonresponse。

即使受訪者接受訪問，也可能拒絕回答某些問題，稱為item nonresponse。

我們無法得知拒絕受訪者的基本資料，所以只能盡量避免受訪者拒答，例如找比較熟練的訪員接觸受訪者，或者提供誘因請受訪者回答問卷，但是因為民眾越來越重視隱私，學術研究也沒有公權力，仍然無法完全避免拒訪。

我們也要避免受訪者拒答，例如題目不能設計的太困難或者太敏感、問卷不能太長等等，或者提供誘因希望受訪者回答全部的問題，但是礙於研究倫理以及個人隱私，還是很難避免拒答。

分析問卷時，我們通常會刪掉拒答該題目的受訪者，無可避免地會產生偏誤，得到的估計的不確定性會提高，所以最好有一定的樣本數，萬一必須刪掉拒答的受訪者，也還能分析。

2 描述統計

以最有效率的方式描述量化或類別變數的重要特徵。
分析單位可以是個體或是群體。
個體的描述統計，例如：個人的性別、教育程度、滿意度、收入、政治態度、血糖指數\(\ldots\)等等。
群體的描述統計，例如：各國的人均所得、房價中位數、收入不均指數（吉尼係數）、生育率\(\ldots\)、平均每千人的醫生人數等等。
個體的集合可以構成集體，例如個人所得的平均是人均所得，個別房價的中位數可以用來代表一個地區的房價水準等等。但是出生率是該年出生人數／該年人口總數\(\times 1000‰\)，不是個體統計而來。
描述統計可分為中央趨勢以及離散程度兩個面向。
可以用表格或者是圖形表示描述統計。

2.1 資料性質

不連續變數或是類別變數
- 全台灣的汽機車數
- 考試分數
- 擲骰子的點數
- 每一個星座的人數
- 臉書上線的人數
- 智商
連續變數
- 溫度
- 高度、重量
- 時間
- 失業率
- 吉尼係數
- 加權指數
如果是不連續變數，用PMF(Probability Mass Function)代表變數中每一個事件的發生機率：

\[ p(x)=P(X=x)\\ 0\le p(x)\le 1 \\ \sum p(x) =1 \]

例如擲一顆公平的骰子，點數的PMF可以表示如下：

\[ p(1)=P(X=1)= \frac{1}{6}\\ p(2)=P(X=2)= \frac{1}{6}\\ p(3)=P(X=3)= \frac{1}{6}\\ p(4)=P(X=4)= \frac{1}{6}\\ p(5)=P(X=5)= \frac{1}{6}\\ p(6)=P(X=6)= \frac{1}{6} \]

如果是連續變數，用PDF(Probability Density Function)描述變數的某一段範圍的發生機率。
假設成人的身高成標準常態分佈，我們可以計算平均值左右一個Z值之間的身高的發生機率如圖 2.1。我們會在後面討論Z值。

# Load ggplot2 package
library(ggplot2)

# Define the range for x values
x_range <- seq(-3, 3, by = 0.01)

# Define the standard normal distribution function
y <- dnorm(x_range, mean = 0, sd = 1)

# Create the data frame for plotting
data <- data.frame(x = x_range, y = y)

# Plot the normal distribution curve
p <- ggplot(data, aes(x = x, y = y)) +
  geom_line() + # Draw the curve
  geom_area(data = subset(data, x >= -1 & x <= 1), aes(x = x, y = y), fill = "blue", alpha = 0.5) +
  labs(title = "Normal Distribution with Shaded Area Between z = -1 and z = 1",
       x = "Z value",
       y = "Density") +
  theme_minimal()

# Display the plot
print(p)

Figure 2.1: 身高常態分佈

而常態分佈的函數可寫成如下：

\[ f(x | \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}} e^{ -\frac{1}{2} \left( \frac{x - \mu}{\sigma} \right)^2 } \]

其中：

\(x\) 代表隨機變數的值
\(\mu\) 代表平均值
\(\sigma\) 代表離散程度
\(\pi\) 則是圓周率(3.1416…)，有常態化(normalize)的作用，\(\frac{1}{\sigma\sqrt{\pi}}\)把資料從\(-\infty\)到\(\infty\)的在曲線底下，而且面積等於1。
\(e\) 自然對數的底，也是一個無理數，約等於2.71。\(e\)可以用R表示為：

e <- exp(1)
e
## [1] 2.718

如果我們把e的平方、次方、立方…連起來，用圖2.2表示：

#simulate
n <- 10
e_powers <- exp(-n: n)
#data
df <- data.frame(x = c(-n:n), y = e_powers) 
#graphic
xg <- ggplot2::ggplot(data = df, aes(x = x, y = y)) +
  geom_line(col = '#CC11EE') +
  geom_point(col = '#1100FF')
xg

Figure 2.2: 自然對數的底

可以看出\(e\)有緩慢上升的趨勢，與常態分佈有關。
假設政大學生身高成常態分佈，平均值是171，標準差是7公分。今天抽樣統計，得知樣本當中身高最矮是160，最高175，請問這套樣本的身高在全體政大學生的當中，出現機率有多高？

# Probability that a randomly selected adult male is between 168 cm and 182 cm
lower_bound <- 160
upper_bound <- 175

# Mean and standard deviation
mean_height <- 171
std_dev <- 7

# Calculating the probabilities
prob_lower = pnorm(lower_bound, mean = mean_height, sd = std_dev)
prob_upper = pnorm(upper_bound, mean = mean_height, sd = std_dev)

# Probability of a randomly selected adult male being between 168 cm and 182 cm
prob_between = prob_upper - prob_lower
prob_between
## [1] 0.6581

0.65代表65%的機率，也就是說這套樣本約代表政大65%的學生。圖2.3顯示樣本落在母體的中間偏左。
如果換成政大學生智商的分佈，我們可以假設成常態分佈嗎？

# Mean and standard deviation
mean_height <- 171
std_dev <- 7
lbound <- 160
ubound <- 175


# Generate data
height <- rnorm(1000, 171, 7)

# Create a data frame with the density values
density <- data.frame(x = seq(min(height), max(height), length.out = 1000),
                      y = dnorm(seq(min(height), max(height), length.out = 1000), mean(height), sd(height)))

# Plot the density curve
hp <- ggplot(density, aes(x = x, y = y)) +
  geom_line() +
  geom_ribbon(data = subset(density, x >= 160 & x <= 175),
              aes(ymax = y, ymin = 0), fill = "#C11100", alpha = 0.3) +
  labs(title = "Normal Distribution with Shaded Area",
       x = "Height", y = "Density")
hp

Figure 2.3: 身高分佈

2.2 描述統計方式

例如UsingR::alltime.movies有電影的票房資料。我們可以找出這筆資料的中央趨勢，例如計算2000年（含）前後的票房平均數：

 movies<-UsingR::alltime.movies; attach(movies)
 dt1 <- movies[which(Release.Year < 2000),]
 dt2 <- movies[which(Release.Year >= 2000),]
 cat("mean before 2000:", mean(dt1$Gross), "\n")
## mean before 2000: 243.5
 cat("mean after 2000:", mean(dt2$Gross))
## mean after 2000: 234.6

另外一個分析方式為：

movies <- movies %>% mutate (D = ifelse (Release.Year < 2000, 'Before 2000', 'After 2000')) %>% 
  mutate (D = as.factor(D)) %>% 
  mutate (D = factor(D, levels=c("Before 2000", "After 2000"))) %>% 
  dplyr::select (D, Gross)
tblmovies<- movies %>% 
  tbl_summary(
            by = D, 
              type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean}",
      "{median} ({p25}, {p75})",
      "{min}, {max}"
    ),
    missing = "no",
        digits = all_continuous() ~ 2
  ) %>% 
  modify_header(label ~ "**Release Year**") %>% 
   modify_caption("**Boxoffice and Release Year**") %>%
  bold_labels()
  
tblmovies

Table 2.1: **Boxoffice and Release Year**
Release Year	Before 2000, N = 50	After 2000, N = 29
Gross
Mean	243.46	234.55
Median (IQR)	216.50 (181.75, 257.75)	215.00 (190.00, 260.00)
Range	172.00, 601.00	176.00, 404.00

圖 2.4顯示，2000年的票房集中在左邊，2000年後有部分電影票房出現在右邊。

library(sjPlot)
H.movies<- movies %>% 
  group_by(D) %>% 
  plot_frq(Gross, type = "histogram", show.mean = TRUE, normal.curve = TRUE) %>% 
  plot_grid()

Figure 2.4: 2000年前後的電影票房分佈

H.movies

TableGrob (2 x 1) “arrange”: 2 grobs z cells name grob 1 1 (1-1,1-1) arrange gtable[layout] 2 2 (2-2,1-1) arrange gtable[layout]

如果是不連續變數，我們改用長條圖，例如在dplyr有一筆starwars的資料，我們想知道電影裡面角色的性別：

p.sex <- starwars %>% sjPlot::plot_frq(sex) 
p.sex

Figure 2.5: Eye Color of Characters in Star Wars

可以看出最多的性別是男性，其次是女性。括號內的百分比代表相對的比例，也就是：

\[ p_{i}=\frac{freq_{i}}{N}\times\frac{100}{100} \]

男性角色佔了全部角色的72.3個百分比。在解讀兩個類別的差距時，要說差距多少百分點。例如男性角色比女性角色多了53個百分點。

如果用沒有畫圖的語法來分析，可以用janitor::tabyl這個函數：

t.sex <- starwars %>% janitor::tabyl(sex) 
t.sex
##             sex  n percent valid_percent
##          female 16 0.18391       0.19277
##  hermaphroditic  1 0.01149       0.01205
##            male 60 0.68966       0.72289
##            none  6 0.06897       0.07229
##            <NA>  4 0.04598            NA

或者用\(\texttt{prop.table}\)函數表示次數分佈：

table(starwars$sex) 
## 
##         female hermaphroditic           male           none 
##             16              1             60              6
prop.table(table(starwars$sex))
## 
##         female hermaphroditic           male           none 
##        0.19277        0.01205        0.72289        0.07229

電影中的男性角色佔了7成。

如果變數的類別太多，我們可以改成橫的長條圖，同樣可以觀察最多次數的類別如圖2.6：

library(ggplot2); library(dplyr)
eyecolors <- unique(starwars$eye_color)
p.eye <- starwars %>% ggplot2::ggplot(aes(x = eye_color, fill = eye_color)) + 
   geom_bar() +
  scale_fill_manual(values= c('black','blue','#AAAAFF', 'brown',
                              'gray20', 'gold', 'greenyellow','#BB2211',
                              'orange','pink', 'red', '#CCEEFF', '#EE11EE',
                              'white', 'yellow')) +
  ggplot2::labs(caption = "Star Wars",   y = 'Frequency', x = 'Eye Color')+
  theme(legend.position = 'none') +
  coord_flip() 
p.eye

Figure 2.6: 星際大戰角色的眼睛顏色

如何比較同一套樣本中不同類別的次數？例如1,000個選民當中，其中有350人支持興建核四，有390人反對，其他人無意見，請問反對的比例顯著高於支持的比例嗎？

我們可以先忽略無意見的人，也就N=350+390=740。然後我們計算390/740=0.527是否顯著大於0.5，如果是的話，那就是支持者過半數，也就不需要比較反對多還是支持多。如果不是，那麼我們比較支持者的比例是否高於反對者比例。也就是說，同一套計算方法，但是比較不同的比例。

n <- 390+350
p1 <- 390/n
p=0.5
s1 <- sqrt(p*(1-p)/n)
z <- (p1-p)/s1; z

[1] 1.47

1-pnorm(z)

[1] 0.07072

p2 <-350/n
s2 <- sqrt(p2*(1-p2)/n)
z2 <- (p1-p2)/s2; z2

[1] 2.945

因為z1<1.96，表示p1並沒有顯著大於0.5，所以支持者並沒有過半。但是z2>1.95，所以支持者明顯比反對者來得多。

如果一個民調題目有超過兩個選擇，每個選擇的類別應該是多元名目分配(multinomial distribution)，兩個類別的比例差距的抽樣誤差應該符合多元名目分配的中央極限定理，得到的抽樣誤差為以下這個公式:

\[ 1.96*\sqrt{\frac{\hat{p_{1}}+\hat{p_{2}}-(\hat{p_{1}}+\hat{p_{2}})^2}{N}} \]

當樣本數是1,000時，除非兩個類別的差距非常大，不然95%信心水準的抽樣誤差大概是5%，樣本數要到3000，抽樣誤差才會降到3%。

p1<-0.36; p2<-0.40; n=1000;
s.error = 1.96*sqrt((p1+p2-(p1-p2)^2)/n); s.error

[1] 0.05398

有四種中央趨勢的統計：

眾數
中位數
百分位數
平均數
從這四個統計值，可以大致判斷資料的分佈狀況。

2.3 中央趨勢：眾數(mode)

眾數適用於間斷變數，例如性別、地區、族群等等，不適用於連續變數。
眾數的定義是發生最多次的那一個類別，例如哪一個節目最多人看。眾數有可能超過一個。
相對於其他類別，眾數所在的類別可以代表較可能發生的事件。如果知道眾數所在的類別，可以用這個類別去猜測或是代表資料以外的事件。
例如已知多數的警察是男性，我們如果隨機抽出一位警察，應該會猜該受訪者是男性。但是我們仍然有 \(100\times (1-m)\%\) 的機會犯錯，\(m\) 代表已知警察為男性的比例，\(1>m>0\)。如果有其他資訊，我們可以降低 \(100\times (1-m)\%\)。
例如已知某一個國家的小學教師之中有6成是女性，\(m=0.6\)，我們有\(100\times (1-0.6)\%=40\%\)誤認某一位老師是女性的可能。
R的mode()函數會回傳向量儲藏資料的性質，並不會告訴我們眾數。例如我們讀了一筆以SPSS格式儲存的民調資料：

b2<-here::here('data','PP0797B2.sav')
dt <- haven::read_spss(b2)
mode(dt$Q1)
## [1] "numeric"
table(dt$Q1)
## 
##   1   2   3   4  95  96  97  98 
## 617 684 443  91  10  57  52 104

我們可以自己寫一個函數來得到眾數，首先我們創造一個向量，呈現變數的表格，然後用names()找出這個表格的首行，進一步篩選首行的元素，條件為該表格的最大值，符合這個條件的就是該變數的眾數：

 tmp <- table(as.vector(dt$Q1))
 tmp
## 
##   1   2   3   4  95  96  97  98 
## 617 684 443  91  10  57  52 104
 names(tmp)
## [1] "1"  "2"  "3"  "4"  "95" "96" "97" "98"
 names(tmp)[tmp == max(tmp)]
## [1] "2"

結果是第二類有最多的數，所以眾數等於2。

☛請練習用ISLR套件中的Carseats資料，找出US變數的眾數所在的類別。

2.4 中央趨勢：中位數(Median, Md, M)

在一個依序排列的數列中，位於中央的數稱為中位數。50\(\%\)的數比中位數大，50\(\%\)的數比它小。通常中位數以\(M\)表示。
中位數與另一個分佈的統計有關：百分位數(percentile)，中位數是第50個百分位數(50th percentile)。
由於中位數本質是排序，與數與數之間的距離無關，所以中位數不會受到極端數值的影響，比較能反映數列的中心位置。但是中位數不適合代數的演算。
中位數可用來表示收入、房屋年齡、房屋坪數、房價，例如2016年我國工業及服務業每人每月業薪資中位數為4萬853元，2019年工業及服務業受僱員工全年總薪資（含經常性與獎金等非經常性薪資）中位數則為49.8萬元（平均每月約4.2萬元），較2018年增加1.64%（資料來源：中華民國統計資訊網）。。
例如：內政部營建署調查公布的房價負擔能力指標，包含「房價所得比」與「貸款負擔率」兩項，「房價所得比」計算公式為「中位數住宅總價÷家戶年可支配所得中位數」，代表需花多少年的可支配所得才買到一戶中位數住宅總價，數值越高表示房價負擔能力越低。
例如在UsingR的套件中，Boston這筆資料有房價中位數的變數 medv，我們用散佈圖表示房價中位數與生師比ptratio以及低社會地位人口比例lstat的關係。圖 2.7顯示，生師比越低、低社會地位人口比例越低，房價的中位數越高。

library(ggplot2)
ggplot(data=MASS::Boston, aes(y=medv, x=ptratio)) +
       geom_point(aes(color=lstat))

Figure 2.7: 波士頓各區的房價中位數與生師比及低社會地位人口比例散佈圖

2.4.1 中位數計算方式

當數列的數目是奇數，中位數是第\(\frac{n+1}{2}\)的數。
如果個數是偶數的資料數列，中位數是\(\frac{a+b}{2}\)，\(a\)、\(b\)是第\(\frac{n+1}{2}\)的數相鄰的兩個數。
例如：1到10，中位數是5.5。0到10，中位數是5。

A<-c(1:10); B<-c(0:10)
median(A); median(B)

[1] 5.5 [1] 5

☛請問studentsfull.txt這筆資料中，學生的中位數成績是多少？

2.4.2 從中位數與平均數看資料

我們用常態分佈模擬資料然後畫圖2.8，發現平均數在中位數的左邊：

# Generate some random data
set.seed(02138)
data <- c(rnorm(500, mean = 40, sd = 12), rnorm(500, mean=60, sd=5))

# Calculate median, mean, and mode
median_val <- median(data)
mean_val <- mean(data)

# Density
plot(density(data), col = "black",  
             main = "Distribution with Central Tendency Measures", 
             xlab = "Value")

# Add lines for median, mean, and mode
abline(v = median_val, col = "red", lwd = 2, lty = 2)
abline(v = mean_val, col = "blue", lwd = 2, lty = 2)


# Add legend
legend("topright", legend = c("Median", "Mean"), col = c("red",  "blue"), lwd = 2, lty = 2)

Figure 2.8: 平均數在中位數的左邊

可以看到分佈向右偏，一半的數在平緩的左邊，一半的數在高聳的右邊。平均數在左邊，代表左邊的數雖然分散，但是總和等於右邊的總和，
我們用均等分佈模擬資料，畫圖2.9表示平均數在中位數的右邊：

# Generate some random data
set.seed(02138)
data <- c(runif(600, 10, 40), runif(400, 65, 250))

# Calculate median, mean, and mode
median_val <- median(data)
mean_val <- mean(data)

# Density
plot(density(data), col = "black",  
             main = "Distribution with Central Tendency Measures", 
             xlab = "Value")

# Add lines for median, mean, and mode
abline(v = median_val, col = "red", lwd = 2, lty = 2)
abline(v = mean_val, col = "blue", lwd = 2, lty = 2)


# Add legend
legend("topright", legend = c("Median", "Mean"), col = c("red",  "blue"), lwd = 2, lty = 2)

Figure 2.9: 平均數在中位數的右邊

跟圖2.9相反，分佈向左偏，平均數在中位數的右邊。

2.5 中央趨勢：四分位數(quantile)

四分位數是數列分成四份之後的三個點：25分位、50分數、75分為其中的25與75分位數。
對於數列的分佈有不同的假設，就有不同計算百分位數的方式。
四分位數是依序排列觀察值，分成四等份的分位數\(Q_{i}\)，\(i=\{1,2,3\}\)，\(Q_{1}\)代表有\(\frac{1}{4}\)的觀察值小於\(Q_{1}\)，\(Q_{3}\)代表有\(\frac{3}{4}\)的觀察值小於\(Q_{3}\)。
例如資料為：\(X=(1, 1001, 1002, 1003)\)
25 百分位所在位置\(=\frac{4\times 25}{100}=1\)。因此 25百分位為 1。
50 百分位所在位置為：\(\frac{4\times 50}{100}=2\)。因此 50百分位為 1001。
75 百分位所在位置為：\(\frac{4\times 75}{100}=3\)。因此 75百分位為 1002。

X <- c(1, 1001, 1002, 1003)
qd <- quantile(X, c(.25, .5, .75), type=1)
knitr::kable(qd, col.names = c('quantiles', 'value')) %>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	1
50%	1001
75%	1002

例如：11位大學生的手機資費如下：195,220, 250,250,305,311,350,371,420,473,650。
\(Q_{1}=\frac{11}{4}=2.75\)。進位之後取第3個數，得到250。
\(Q_{2}=\frac{2\times11}{4}=5.5\)。進位之後取第6個數，得到311。
\(Q_{3}=\frac{3\times11}{4}=8.25\)。進位之後取第9個數，得到420。
R提供9種計算方法，前面3種適用於間斷變數，後面6種則是連續變數。我們以第1種方法計算。

m <- c(195,220, 250,250,305,311,350,371,420,473,650)
qd <- quantile(m, c(0.25, 0.5, 0.75), type = 1)
knitr::kable(qd, col.names = c('quantiles', 'value'))%>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	250
50%	311
75%	420

例如我們模擬兩筆資料，平均數相同，但是離散程度不同，以同樣的方法計算25, 50, 75分位，會得到不同的25及75分位。這個例子顯示即使中位數相同，但是資料分佈不同，25及75分位就不同。

set.seed(1000)
tmp1<-rnorm(100, 40, 10)
qmp1 <- quantile(as.integer(tmp1), c(0.25, 0.5, 0.75), type=1)
knitr::kable(qmp1, col.names = c('quantiles', 'value'))%>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	34
50%	40
75%	45

set.seed(1000)
tmp2<-rnorm(100, 40, 14)
qmp2 <- quantile(as.integer(tmp2), c(0.25, 0.5, 0.75), type=1)
knitr::kable(qmp2, col.names = c('quantiles', 'value'))%>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	31
50%	40
75%	47

2.6 中央趨勢：百分位數(percentile)

把資料由小排到大，第\(i\)個百分位數表示(100-\(i\))%的數比它大，\(i\%\)的數比它小。可以表示資料的集中與分散。也被稱為百分等級(percentile rank)。例如pr99是286分。
利用累積相對次數，用1％, 2％, 3％,\(\ldots\), 99％將資料均分成100等份，中間99個分割點所得到對應的數值，稱為該資料的第1、2、3…、99百分位數。
可以是實際存在的數，也可以是計算所得。
有數種計算方式，應該根據資料的分佈（或假設）而選擇計算方式。其中一種百分位數的計算公式為：

\[m_{i}=100\cdot \frac{i}{n}\]

上面的公式中，\(m\)變數的\(i\)百分位數等於\(i\)除以\(m\)變數的觀察值總數\(n\)再乘以100。如果\(m_{i}\)不是整數，則\(k\)為該百分位數，且\(m_{i+1}\ge k\ge m_{i}\)。
換句話說，當\(m_{i}\)不是整數，我們可以將\(m_{i}\)無條件進位加1的數當做\(m_{i}\)。
另一種算法是當\(m_{i}\)是整數，則排在第\(m\)位與\(m＋1\)位資料值的算術平均數就是這群資料的第\(k\)百分位數。
我們也可以用另一筆資料驗證手算以及R的結果：

full<-here::here('data','studentsfull.txt')
dt <- read.csv(full,sep="",header=T)
dt$Score<-sort(dt$Score)
dt$Score

[1] 60 62 66 66 69 70 75 77 78 80 80 81 82 83 85 85 86 87 88 88 88 89 91 92 92 [26] 93

dt$Score[floor(length(dt$Score)*0.25)+1]

[1] 75

dt$Score[floor(length(dt$Score)*0.75)+1]

[1] 88

dt$Score[floor(length(dt$Score)*0.9)+1]

[1] 92

qd1 <- quantile(dt$Score, c(0.25, .75, 0.9), type=1)
qd1 <- knitr::kable(qd1, col.names = c('quantiles', 'value'))%>%
  kable_styling(latex_options = "scale_down")
qd2 <- quantile(dt$Score, c(0.25, .75, 0.9), type=4)
knitr::kable(qd2, col.names = c('quantiles', 'value'))%>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	72.5
75%	88.0
90%	91.4

2.6.1 比較SPSS與R的輸出

R的輸出跟SPSS類似，我們可以加以對照（圖2.10）。SPSS的統計值等於是R的quantile()的第六種計算方式。
例如有一筆34位學生的成績資料，我們計算25, 50, 75, 90百分位的數字：

scores<-c(15, 22, 26, 32, 33,36, 36, 41, 42, 44,
44, 45, 47, 48, 61,63, 63, 65, 65, 65,
66, 66, 68, 69, 70,71, 74, 74, 76, 77,
78, 78, 80, 85)
qscore <- quantile(scores, c(.25,.5,.75,.9), type=6)
knitr::kable(qscore, col.names = c('quantiles', 'value')) %>%
  kable_styling(latex_options = "scale_down")

quantiles	value
25%	41.75
50%	64.00
75%	71.75
90%	78.00

v1<-here::here('Fig','v1_quantile.png')
knitr::include_graphics(v1)

Figure 2.10: 四分位統計

2.7 中央趨勢：平均數

平均數衡量資料的中心位置，可以想成是觀察值的平衡點：比平均值大的數的總和等於比平均值小的數的總和的絕對值。
平均數會受到極端值的影響，可以用trim刪除若干百分比的數。
可分為算術平均數跟加權平均數。

2.7.1 算數平均數：

算術平均數(arithmetic mean)的公式如下：\[\bar{y}=\frac{\sum y_{i}}{n}\]
例如： \[y={6, 7, 8, 8, 9, 10, 13, 15, 16, 45}\]
平均數為： \[\bar{y}=\frac{\sum (6+7+\cdots , +45)}{10}=13.7\]
如果\(x_{i}=y_{i}+10\)，請問平均數\(\bar{x}\)會比\(\bar{y}\)大、小、不變？
如果\(x_{i}=10\times y_{i}\)，請問平均數\(\bar{x}\)會比\(\bar{y}\)大、小、不變？

2.8 加權平均數

在不知道個別觀察值，只知道分組的個案數跟平均數，我們可以假設觀察值分為\(k=1\cdots k\)個組，每一組有\(y_{1}\),\(y_{2}\),\(\ldots\) 人，每一組平均數為\[\bar{y_{1}}, \bar{y_{2}},\cdots\], 則全體的平均數為：

\[\bar{y}=\frac{\sum n_{k}\cdot \bar{y_{k}}}{n}\]

換句話說，每一組的平均數乘以每一組在全部資料中佔的比例，類似加權，得到的平均數和就是全體樣本的平均數。
例如有三個空氣品質的觀測站的資料，要計算總平均數，首先從總和除以全部個案數計算：

A<-list(station1=c(25, 33, 44),
        station2=c(43, 66, 78, 81),
        station3=c(90, 76, 105, 110, 121))
#n of each group
group.n=sapply(A, length); group.n
## station1 station2 station3 
##        3        4        5
#n of data
total.n=sum(sapply(A, length));
#proportion of each group
group_p <- group.n/total.n
#mean of each group
submean=sapply(A, mean); cat("air pollution of each station:", submean,"\n")
## air pollution of each station: 34 67 100.4
cat("Average air pollution=", sum(group_p*submean),"\n")
## Average air pollution= 72.67

#sum of data
totalair=sum(sapply(A, sum));
cat("Sum of air pollution=",totalair,"\n")
## Sum of air pollution= 872
cat("average air pollution=", totalair/total.n)
## average air pollution= 72.67

sapply()函數可套用函數在列表的每一個向量。

2.9 偏態(skewness)

偏態表示變數的分佈的對稱程度。圖2.11顯示正、負偏態的型態。

有偏態時須注意平均值是否會誤導。
正偏：右邊的尾巴較左邊長，眾數<中位數<平均數，偏態係數大於0。
負偏：左邊的尾巴較右邊長，眾數>中位數>平均數，偏態係數小於0。不過眾數的位置不一定在最右邊。
常態分佈的偏態值＝0
樣本偏態值\(=\frac{n}{(n-1)(n-2)} \frac{\sum (x_{i}-\bar{x})^3}{s^3}\)
\(s\)是標準差\(s=\sqrt{\frac{\sum (x-\bar{x})^2}{n}}\)

spss<-here::here('Fig','week3_skewness.jpg')
knitr::include_graphics(spss)

Figure 2.11: 偏態統計

2.9.1 R與Stata以及SPSS的比較

偏態有多種計算方式，R的計算公式1與2分別與Stata以及SPSS得到的結果相同。
以學生的寫作成績為例，偏態分別是-0.47以及-0.48，表示平均的寫作成績小於中位數：

library(foreign)
hs<-here::here('data','hsb2.dta')
hsb2<-read.dta(hs)
library(e1071)
cat('skewness of writing:', skewness(hsb2$write, type=1), '\n')

skewness of writing: -0.4784

cat('skewness of writing:', skewness(hsb2$write, type=2))

skewness of writing: -0.482

首先是Stata偏態計算結果，如圖2.12：

stata<-here::here('Fig','write_stata.png')
knitr::include_graphics(stata)

Figure 2.12: Stata的偏態統計

再來是SPSS偏態計算結果，如圖2.13：

spss<-here::here('Fig','write_spss.png')
knitr::include_graphics(spss)

Figure 2.13: SPSS的偏態統計

2.9.2 模擬資料

假設模擬三組資料，計算偏態如下：

# Sample data
set.seed(02138)
A <- c(runif(700, 20, 30), runif(300, 30, 70))
B <- rnorm(1000,  30, 4)
C <- c(runif(200, 0, 80), rnorm(800, 50, 10))
cat('A skewness:', skewness(A, type=1), '\n')

A skewness: 1.396

cat('B skewness:', skewness(B, type=1), '\n')

B skewness: 0.01334

cat('C skewness:', skewness(C, type=1), '\n')

C skewness: -0.7783

計算結果顯示：A 是正偏，B是常態分佈，沒有偏態，C是負偏。
用表??呈現三筆資料的平均值與中位數，然後計算眾數。A資料的眾數最小，B資料跟平均數、中位數相近，B資料的眾數也接近平均數、中位數

library(gtsummary)
library(dplyr)
library(e1071)
# Create the dataset
tmp <- data.frame(A = A, B = B, C = C)

# Create tbl_summary object
tblm1 <- tmp %>% 
  tbl_summary(
              type = all_continuous() ~ "continuous2",
    statistic = all_continuous() ~ c(
      "{mean}",
      "{median}"
    ),
    missing = "no",
        digits = all_continuous() ~ 2
  ) %>% 
  bold_labels() 
tblm1

Characteristic	N = 1,000
A
Mean	32.30
Median	27.04
B
Mean	29.96
Median	29.87
C
Mean	48.15
Median	49.55


#mode
varmode <-  function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

cat('Mode of A:', varmode(A), '\n')

Mode of A: 25.12

cat('Mode of B:', varmode(B), '\n')

Mode of B: 26.72

cat('Mode of C:', varmode(C), '\n')

Mode of C: 50.36

#tblm2 <- tmp %>% 
 # tbl_custom_summary(
 #stat_fns = ~varmode,
 #statistic = ~'{mode}'
  #  )
  

#tbl_merge_ex2 <-
#  tbl_merge(tbls = list(tblm1, tblm2)) %>%
#  modify_spanning_header(everything() ~ NA_character_)

我們可以用橫的盒型圖如圖2.14表示分佈：
- A 資料之中，有一半（25分位到75分位）集中在20到40之間，中位數大概接近30，分佈右邊有很多觀察值。
- C 資料的落在右邊，但是左邊也有一些分散的觀察值。B資料則集中在30左右。

data <- data.frame(group = rep(c("A", "B", "C"), each = 1000),
             value = c(A, B, C))
# Define colors for each group
colors <- c("#EE22AA", "#EBC311", "#C1F210")

# Draw horizontal box plot
data$group <- ordered(as.factor(data$group), levels = c('C','B','A'))
boxplot(value ~ group, data = data, horizontal = TRUE, col=colors,
        ylim=c(0, 80))
# Adjust y-axis ticks
axis(2, at = seq(0, 80, by = 10))

Figure 2.14: 表示偏態的箱型圖

畫密度圖2.15對照：

x <- seq(min(A), max(A), length = 1000)
df <- data.frame(x=c(x, x, x), y=c(sort(A), sort(B), sort(C)), 
   dataset = c(rep('A', 1000), rep('B', 1000), rep('C', 1000)))
ggplot2::ggplot(df, aes(x = y, group = dataset)) +
    geom_density(aes(fill = dataset), alpha = 0.5) +
    scale_fill_manual(values=c("#C1F210", "#EBC311", "#EE22AA"))

Figure 2.15: 3個模擬資料的密度圖

2.10 峰度 (kurtosis)

峰度的用途是測量資料的分佈是高聳或是平坦的程度。
越平坦則兩邊尾部越長，越高聳則是靠近平均值的部份越集中。
計算峰度的公式為： \[\frac{m^4}{m^2_{2}}-3\] \[s_{2}=\sum (x_{i}-\bar{x})^2\] \[s_{4}=\sum (x_{i}-\bar{x})^4\] \[m_{2}=\frac{s_{2}^2}{n}\] \[m4=\frac{s_{4}}{n}\]
不同統計軟體計算峰度的公式略有不同，如果用R，可以用e1071這個套件裡面的kurtosis，用type指令選擇2可得到跟SPSS一樣的答案。
Stata的計算峰度公式為 \[\frac{m^4}{m^2_{2}}\]
樣本數目越大，理論上各種計算方式的結果越接近。

2.10.1 使用語法計算峰度

#m4
s4=sum((hsb2$write-mean(hsb2$write))^4)
m4=s4/200; print(m4)
## [1] 17889
#m2
s2=sum((hsb2$write-mean(hsb2$write))^2)
m2=s2/200; print(m2)
## [1] 89.39
#kurtosis
cat('kurtosis=', m4/m2^2)
## kurtosis= 2.239

2.10.2 比較R、Stata與SPSS

Stata的計算方式是\(\frac{m^4}{m^2_{2}}\)，而SPSS則是\(\frac{m^4}{m^2_{2}}-3\)。
首先是兩種R計算結果：

kurtosis(hsb2$write, type=1)+3

[1] 2.239

kurtosis(hsb2$write, type=2)+3

[1] 2.25

其次是Stata，如圖2.16：

stata<-here::here('Fig','write_stata.png')
knitr::include_graphics(stata)

Figure 2.16: Stata的峰度統計

最後是SPSS，如圖2.17：

spss<-here::here('Fig','write_spss.png')
knitr::include_graphics(spss)

Figure 2.17: SPSS的峰度統計

以上的結果顯示我們可以互相比較三種軟體計算的結果，R的彈性比較大。

2.11 峰度範例

模擬常態分佈與t分佈的資料，計算常態分佈的峰度為0，t分佈則是74。

# Generate random samples from the normal distribution
normal_dist <- rnorm(1000)
e1071::kurtosis(normal_dist, type=1)
## [1] 0.02436
# Generate random samples from the t-distribution with 3 degrees of freedom
t_dist <- rt(1000, df = 3)
e1071::kurtosis(t_dist, type=1)
## [1] 46.38

用長條圖2.18表示t分佈：

hist(t_dist, freq = FALSE, main = "t-Distribution (df = 3)", 
     xlab = "Value", ylab = "Density", col = "lightgreen", ylim=c(0, 0.33),
     breaks=40)
lines(density(t_dist), col = "#AA11AA", lwd = 2)

Figure 2.18: T分佈模擬PDF

3 離散程度

離散指的是資料點與中央趨勢之間的差異，包括中位數以及平均數。
資料離散有可能本身就很離散，或者是因為測量誤差。
抽樣誤差可能導致我們觀察到離散程度大的資料。
抽樣時要小心離散是否與特定群體有關，例如金融業可能因為工作性質是理財顧問、銀行行員或者基金經理人而有很大的收入差異。
以下依照範圍、四分位差、變異數、標準差依序介紹。

3.1 範圍

範圍(range)：最大值及最小值的差距，也稱做全距。
若是常態分佈，範圍約等於六個標準差。

range(hsb2$write)

[1] 31 67

range(hsb2$math)

[1] 33 75

平均數相同，範圍可能不同；有的變數比較離散。
全距相同，但是離散程度可能不同，請見圖 3.1。

par(bg="#33110011")
x <- seq(-3, 3, length=100)
hx <- dnorm(x)
plot(x, hx, type="l", lty=2, xlab="x value", lwd=2,
  ylab="Density", main="Comparison of t Distributions", ylim=c(0,0.6))

y <- seq(-3, 3, length=600)
curve(dnorm(x, 0, 0.75), type="l", add=T, lwd=2, col="red")

Figure 3.1: 相同全距離散程度不同的機率密度圖

全距容易受到最大值與最小值的影響；最大值會拉大全距，離散程度變大。

3.2 四分位距(IQR)

第一個四分位跟第三個四分位之間的差距。不受到極端值的影響。
可用在箱型圖或盒鬚圖，例如模擬一筆資料，然後畫箱形圖如 3.2：

library(ggplot2)
sv<-c(15, 22, 26, 32, 33,36, 36, 41, 42, 44,
44, 45, 47, 48, 61,63, 63, 65, 65, 65,
66, 66, 68, 69, 70,71, 74, 74, 76, 77,
78, 78, 80, 85)
quantile(sv, c(.25,.5,.75), type=6)
##   25%   50%   75% 
## 41.75 64.00 71.75
dt <- data.frame(scores=sv)
ggplot(data=dt, aes(y=scores)) +
      geom_boxplot(fill="#FF22EE11")

Figure 3.2: 學生成績盒鬚圖

可以看到，盒鬚圖的下緣約等於\(Q_{1}\)，上緣則是\(Q_{3}\)，中間的線則是中位數。盒鬚圖的中間部分是四分位距，也就是71.75-41.75=30。

☛請計算MASS::Animals這筆資料中的腦容量的四分位距。

3.3 樣本變異數：母體變異數的無偏估計

每一個觀察值與平均數的差距稱為變異數，用來衡量連續變數資料離散的程度。
母體的變異數公式為：

\[\sigma^2=\frac{\sum (X-\mu)^2}{n}\]

樣本變異數代表樣本觀察值與平均數之間的差距，公式為：

\[S^2=\frac{\sum (X-\bar{X})^2}{n-1}\]

標準差的分母是\(n-1\)是為了避免低估變異數。
樣本變異數的平方根為標準差：\(s=\sqrt{S^2}\)
如果樣本來自二元分佈，在\(n\)次的實驗或者樣本當中，\(p\)為事件發生的機率，平均數為\(p\)，變異數為：

\[n\cdot p\cdot(1-p)\]

標準差為：

\[\sqrt{\frac{np(1-p)}{n-1}}\]

3.3.1 範例一

某10家上市公司的股票價格分別為13.5, 22.2, 31.2, 15.2, 20.3, 21.9, 18.3, 25.3, 21.3, 19.8。請問這些公司的股價的變異數是多少？
如果這10家的股票價格都漲了1塊錢，請問變異數變大、變小或者不變？為什麼？
如果這10家的股票價格都漲了1倍，請問變異數變大、變小或者不變？為什麼？

3.3.2 範例二

某戶人家收集過去12個月的水費帳單，想知道每個月水費有沒有超過300元，發現7, 8, 9月之外，其他月份都沒有超過。請問水費有無超過300元的變異數是多少？
如果該住戶發現其實10月份的水費也超過300元，請問變異數變大、變小還是不變？為什麼？
如果該住戶統計到隔年1月，發現還是3個月的水費超過300元，請問變異數變大、變小還是不變？為什麼？

3.4 常態分佈的樣本標準差

如果隨機變數屬於常態分配，大部分的值應該聚集在平均值加減一個標準差的範圍內，因此，樣本標準差的大小特別重要。
當樣本來自於常態分配的母體，利用微積分可求出平均數的加減1個標準差包含約\(68\%\)的樣本。2個標準差包含約\(95\%\)的樣本。3個標準差包含約\(99\%\)的樣本。
先寫程式計算標準差，再用R的sd()函數驗證：

#hsb2 data
file <- here::here('data', 'hsb2.txt')
hsb2 <- read.table(file, header = T)
#variance
v.write<-var(hsb2$write); sqrt(v.write)
## [1] 9.479
#standard deviation
std = function(x) sqrt(var(x))
std(hsb2$write)
## [1] 9.479
#self-defined function
sd<-function(V)sqrt( sum((V - mean(V))^2 /(length(V)-1)))
sd(hsb2$write)
## [1] 9.479

由以上結果可知，一個樣本標準差等於\(\sqrt{\frac{\sum (X-\bar{X})^2}{n-1}}\)。

3.5 標準差的特性

改變樣本的單位，標準差也會改變，例如有一筆容量的資料，除以1000之後從立方公分變成公升，分別求標準差，發現除以1000之後標準差也縮小1000倍了：

H<-c(15000,7000,19000,3000,15000,19000,4000,12000,
       17000,  9000)
h<-c(15,7,19,3,15,19, 4,12,17, 9)
sd(H); sd(h)
## [1] 5963
## [1] 5.963

加減樣本的值會改變平均值，但是不會改變標準差，因為\[\sum_{i=1\sim n} (x_{i}-\bar{x})\]變成\[\sum ((x_{i}+k)-\overline{x+k})=\sum x_{i}+k-\frac{\sum x}{n}-\frac{nk}{n}=\sum (x_{i}-\bar{x})\]

4 標準常態分佈的Z值

我們用一個分數的正負判斷某個樣本的連續變數是否高於平均數，用絕對值判斷距離平均數的差距有多遠，這就是Z分數或者Z值。
Z值的優點是將原始分數轉換為標準分數後，不同單位的數量或者不同測量的分數，均可直接比較，不會受測量困難程度的影響。
如果轉換原始分數轉換為Z值得到小數或負數，也可以乘以一個常數，再加一常數，變成整數。也就是整個分佈向右邊移動。

透過計算\(Z\)值，常態分佈可以轉換為標準常態分佈，。

\[f(Z)=\frac{1}{2\pi}e^{-\frac{Z^2}{2}}\\ Z=\frac{X-\mu}{\sigma}\]

標準常態分佈的性質為：

\[Z\sim N(0,1)\]

4.1 Z值或分數

\(Z\)值可幫我們瞭解觀察值在資料中的相對位置。計算\(Z\)的公式為：

\[Z=\frac{x-\bar{x}}{s}\] - s是標準差。

母體資料的標準化觀察值以比較觀察值與平均值之間的距離則是：

\[Z=\frac{X-\mu}{\sigma}\] \[\sigma\neq 0\]

如果是標準化常態分佈，也就是平均數為0、變異數為1，\(Z\)值大約介於-6到6之間。

pnorm(6, 0, 1)
## [1] 1
pnorm(-6, 0, 1)
## [1] 9.866e-10

\(Z\)值可轉換為百分位，百分位也可轉換為\(Z\)值。例如標準化常態分佈的\(Z\)=-1.96時，面積\(\approx 2.5\%\)。

pnorm(-1.96,0,1)
## [1] 0.025

可畫圖如圖4.1：

curve(dnorm(x),
       xlim = c(-3, 3),
       ylab = "Density",
       #main = "機率密度與區域",
       col='red', lwd=2, xlab='Z')
cord.1x <- c(-3,seq(-3, -1.96,0.01),-1.96)
cord.1y  <- c(0,dnorm(seq(-3, -1.96,0.01)),0)

polygon(cord.1x,cord.1y,col='grey80')

Figure 4.1: 標準常態分佈曲線下的左邊2.5%區域

可以試試看Z=1.96時，落在分佈右邊的2.5%:

curve(dnorm(x),
       xlim = c(-3, 3),
       ylab = "Density",
       #main = "機率密度與區域",
       col='red', lwd=2, xlab='Z')
cord.1x <- c(1.96 ,seq(1.96, 3, 0.01), 3)
cord.1y  <- c(0, dnorm(seq(1.96, 3, 0.01)), 0)

polygon(cord.1x,cord.1y,col='grey80')

Figure 4.2: 標準常態分佈曲線下的右邊2.5%區域

但是\(\texttt{pnorm(1.96, 0, 1)}\)是左邊開始累積的機率，所以1-\(\texttt{pnorm(1.96, 0, 1)}\)才是0.25。

Z = 1.96
1-pnorm(Z, 0, 1)
## [1] 0.025

反過來，\(\texttt{qnorm()}\)顯示若干百分位對應的Z值，例如97.5%對應的Z值是：

q = 0.975
qnorm(q)
## [1] 1.96

又例如：

q = 0.95
qnorm(q)
## [1] 1.645

4.1.1 範例一

★請問在alr4::Heights這筆有關母親跟女兒的身高資料中，請問母親身高mheight介於63英吋（約160公分）與65英吋（約165公分）的比例有多少？

要求\(P(Z\geq z^{*})\)，其中\(z^{*}=\frac{X-\mu}{\sigma}\)，我們先算出\(\mu\)，再計算變異數，得到\(z^{*}\)之後，以pnorm轉換為百分比。

head(alr4::Heights, 4)

mheight dheight 1 59.7 55.1 2 58.2 56.5 3 60.6 56.0 4 60.7 56.8

m.i<-mean(alr4::Heights$mheight)
m.i

[1] 62.45

s.i<-var(alr4::Heights$mheight)
s.i

[1] 5.547

zstar1=(63-m.i)/sqrt(s.i);zstar2=(65-m.i)/sqrt(s.i)
cat("63in=",zstar1,"\n","65in=",zstar2,"\n")

63in= 0.2323 65in= 1.082

pnorm(zstar2, 0, 1)-pnorm(zstar1, 0, 1)

[1] 0.2684

可畫圖4.3表示0.232與1.081之間的區域：

curve(dnorm(x),
       xlim = c(-3, 3),
       ylab = "Density",
       #main = "機率密度與區域",
       col='red', lwd=2, xlab='Z')
cord.1x <- c(0.232,seq(0.232, 1.081,0.01), 1.081)
cord.1y  <- c(0,dnorm(seq(0.232, 1.081,0.01)),0)

polygon(cord.1x,cord.1y,col='grey80')

Figure 4.3: 標準常態分佈曲線下的特定區域

4.1.2 範例二

★ 有一位員工的今年月薪為8.5萬，去年則為8萬。今年的全體薪水標準差為2.3萬，平均值為6.4萬，而去年的全體員工薪水標準差為2萬，平均值為6.2萬。假設員工薪水常態分配。請問該員工月薪相較於全體員工，今年比去年有增加嗎？

求出今年的Z值：

\[z_{1}=\frac{8.5-6.4}{2.3}=0.91\]

求出去年的Z值：

\[z_{2}=\frac{8-6.2}{2}=0.75\] - 因為\(z_{1}\geq z_{2}\)，因此該員工月薪相較於全體員工有增加。

score =78
mu = 67.09
sigma = 11.08
(score-mu)/sigma

[1] 0.9847

5 作業

1. 請計算studentsfull.txt這筆資料中的score中位數、90百分位數以及男性跟女性的平均數：

2. 請計算studentsfull.txt這筆資料中男性跟女性的score平均數以及標準差：

3. 請使用UsingR套件中的faithful資料，請問要看噴泉最少要等幾分鐘？平均要等幾分鐘？最多跟最少等的時間差距多少分鐘？

4. 請用airquality這筆資料的Wind這個變數，計算前後兩個資料點的差異，以分析兩天之間風速的差異。例如：

A <- c(35, 61, 69)
d.A <- c(26, 8)

5. 請問airquality這筆資料的Wind的偏態為何？峰度是多少？

6. 請問在ISLR::College這筆資料中，Private這個變數的樣本標準差是多少？（提示：私立學校設為1，公立學校設為0）

7. 在councilor這筆資料中，請問平均工程預算是多少？樣本標準差多少？

8. 使用2008Election這筆資料，請問馬英九的得票數的25分位數、中位數、75分位數分別是多少？請問25與75分位數之間差別多少？

9. 請打開BES這筆資料，計算表示英國應該離開歐盟(leave)的比例。

10. 同樣的，請回答哪一個投票意向(vote)的次數最多？相對頻率多少？

6 更新講義時間

最後更新時間: 2024-03-22 13:10:45

社會科學統計方法

描述統計

蔡佳泓

3/15/2024