社會科學統計方法
隨機變數、機率原理及描述統計
蔡佳泓
4/1/2022
1 機率
1.1 樣本空間
疫情升三級的機率有多少？明天下雨的機率多高？大樂透中獎的機率有多高？
以大樂透的中獎機率為例，我們可以先計算所有號碼的組合，每一組號碼的中獎機率為： \[ P(A)=\frac{n(A)}{n(S)}=\frac{n(A)}{n} \]
其中A代表頭獎的中獎號碼，n(A)代表中獎號碼的次數，例如有三個人買了同樣的號碼，n(A)=3。n(S)代表所有號碼的組合數，例如大樂透是六個01到49的號碼，六個全對得到頭獎，所以一共有\(49\times 48\times\ldots\times 44=10068347510\)組號碼。
機率也可以定義為無數次實驗之後，某事件A出現的相對次數，\(P(A)=\frac{n(A)}{S}\)。隨機實驗是一個過程，在實驗前我們知道所有可能出現的結果，但是不確定是否會出現，也不知道會出現幾次。一個實驗或是隨機變數中所有可能出現的結果 (outcomes) 所形成的集合稱為樣本空間\(S\)。
例如抽出一張撲克牌，有四種結果：\(S=\{\)紅心, 鑽石, 黑桃, 梅花\(\}\)，總共有52種結果：\(S=\{\)紅心1, 紅心2,, 梅花K\(\}\)。
又例如政大校長選舉第一階段，有三種人可以投同意票：\(S=\{\)學生，教師，行政人員\(\}\)，候選人必須獲得一定比例的同意票才能進入第二階段。
由於每一種可能的結果皆屬於樣本空間，故亦可稱為元素或樣本點。
事件 (event) 為結果所形成的集合，事件為 \(S\) 的子集,以大寫字母 (\(A, B, C, \ldots\)) 表示，可分為以下兩種類型:
1.事件中只有包含一個元素稱作簡單事件 (simple event, 亦稱為樣本點) 例如投擲一顆骰子得點數 3。 2.事件中包含兩個以上的元素稱作複合事件 (compound event)，例如投擲兩顆骰子得點數和為 3 或 5 或 7。
若 A, B 兩事件滿足 \(A\cap B=\emptyset\)時, 則稱此A事件與B事件為互斥事件 (disjoint events)。A, B兩個事件不可能同時發生，這兩個事件中的元素不會互相重疊。
例如學生的年級成為一個樣本空間：\(S=\{\)大一，大二，大三，大四\(\}\)。一名學生不可能同時是大一以及大二的學生，所以屬於大一以及屬於大二是互斥事件。但是一個學生可以同時主修政治學以及會計學，所以主修可能不是互斥事件。
1.2 機率公設
Laplace定義的古典機率公設：
- 對於任意事件 \(A\) 皆滿足：\(0\leq P(A) \leq 1\)。（Non-negativity）
- 樣本空間 \(S\) 的機率等於 1，寫作 \(P(S) = 1\)。（Normalization）
- 如果 \(A_{1}\)，\(A_{2}\)，\(A_{3}\)，為一組有限或者可數的無限的事件且彼此為互斥事件，則這些事件所聯集的機率等於個別事件的機率和等於 \(P(A_{1} \cup A_{2} \cup A_{3} ∪ · · ·) = P(A_{1}) + P(A_{2}) + P(A_{3}) + \cdots\) （Additivity）。
- \(P(\emptyset)=0\)
- 若\(A,B\)為\(S\)中的二事件，且\(A \subset B\)，則 \(P(A) \leq P(B)\)
- 若\(A,B\)為\(S\)中的二事件，則 \(P(A \cup B)=P(A)+P(B)-P(A \cap B)\)
- 若\(A \subset S\)為一事件，則\(P(A')=1-P(A)\)
1.3 實例
- 丟一枚硬幣，樣本空間\(S=\{\)正面，反面\(\}\)，出現正面的事件\(E=\{\)正面\(\}\)
- 丟一個六面的骰子，樣本空間\(S=\{1,2,3,4,5,6\}\)，出現偶數的事件\(A=\{2,4,6\}\)，出現偶數的事件\(B=\{1,3,5\}\)，兩者為互斥事件。
- 丟一枚硬幣兩次，樣本空間\(S=\{HH, HT, TH, TT\}\)，第一次丟擲出現正面的事件為\(E=\{HH, HT\}\)
- 丟一個骰子兩次，樣本空間\(S=\{(i,j):i,j=1,2,\ldots,6\}\)，共有36個元素。兩次骰子點數總和為10為事件之一，也就是\(\{(4,6),(5,5),(6,4)\}\)。
- 抽一張牌，可能是鑽石，也可能是黑桃，抽到這兩者其中之一的機率等於是分別抽到這兩者的機率和，也就是\(P(\text{diamond}\cup \text{spade})=P(\text{diamond})+P(\text{spade})=0.25+0.25=0.5\)
☛請問這兩個函數是否符合機率公設？ \[f(x)=\frac{x-1}{6}\quad x=1,2,3,4\] \[f(x)=\frac{x^2}{2}\quad x=-1,0,1\]
☛提示：(1)\(f(x)\)是否符合\(0\leq f(x)\leq 1\)；(2)代入x到\(f(x)\)，求得\(\sum f(x)=1\)。
1.4 聯合機率 (Joint probability)
- 定義：兩個事件或者兩個隨機變數同時發生的機率，或者是兩個事件交集的機率，可表示為 \(p(A,B)\) 或者是 \(p(A\cap B)\)。如果以隨機變數的質量函數表示，則寫為
\[f(x,y)=P(X=x,Y=y)\]
例如: 一張撲克牌，同時為紅色且為4點的機率為\(P(four, red)=2/52=1/26\)。
美國在2022年出現第一位非裔女性的最高法院大法官：Ketanji Brown Jackson，包括她在內，歷史上116位大法官當中，總共有6位女性大法官，3位黑人大法官，也就是P(Black, Female)=1/116=\(0.8\%\)。
台灣在109年有114,606對結婚，我們統計本國民眾各個年齡層以及性別的結婚人數，發現男性且15至19歲佔全體結婚人數有\(1.6\%\)，也就是P(15_19, Male)=\(1.6\%\)，男性且20至24歲結婚有\(10.6\%\)，依此類推。
file <- '/Users/user/Downloads/opendata109m113.csv'
dt <- read.csv(file, header = T)
library(dplyr); library(janitor)
dt <- dt[-1,]
tab1 <- dt %>% group_by(nation, sex, age) %>%
summarise(n = sum(as.numeric(marry_count))) %>%
filter(nation=='本國') %>%
mutate(total=sum(n)) %>%
mutate (per=n/total)
knitr::kable(tab1,format = 'simple',
booktabs=T,caption="109年本國民眾年齡層以及性別的結婚統計")1.5 獨立事件 (Independent Events)
- 如果事件 \(A\) 發生時與事件 \(B\) 沒有關係，也就是 \(P(A,B)=P(A)\cdot P(B)\), \(P(A|B)\cdot P(B)=P(B)\)，這兩個事件稱為獨立事件。
- 例如：擲一顆骰子兩次，第一次得到6的機率與第二次得到6的機率相互獨立。兩次都得到6的機率是\(\frac{1}{6}\cdot\frac{1}{6}=\frac{1}{36}\)
- 抽一張撲克牌，放回去之後再抽一次，第一次抽到紅心的機率與第二次抽到紅心的機率相互獨立。但是如果抽出不放回，那麼這兩次抽到紅心的機率並不獨立，第二次抽到紅心的機率受到第一次抽到的結果影響。
1.6 互斥事件 (Disjoint Events)
- 兩個事件不可能同時發生，稱為互斥事件，也就是聯合機率為0。兩個互斥事件的機率 \(P(A\cup B)=P(A)+P(B)\)
- 例如，一個學生不可能同時為大一學生，也是大四學生，所以是互斥事件。但是一個學生可以同時為大一學生也是男生。
- 一個生物可能是動物，也可能是植物，但不可能同時是動物跟植物，例如珊瑚是動物，不是植物。
- 例如一個事件是奇數另一個事件是偶數: \(C=\{1,3,5\}, D=\{2,4,6\}\). \(C\) 與 \(D\) 事件為互斥。 \(C\cap B=\emptyset\). \(P(C,D)=0\)
1.7 實例
從52張撲克牌抽到紅色的牌與抽到有臉的牌(J, Q, K)是否為獨立事件?
- 已知抽到紅色牌的機率\(P(Red)=\frac{1}{2}\)
- 已知抽到有臉的牌的機率\(P(Face)=\frac{12}{52}=\frac{3}{13}\)
- 抽到紅色且有臉的牌的機率\(P(R,F)=\frac{6}{52}=\frac{3}{26}\)
- \(\frac{1}{2}\cdot\frac{3}{13}=\frac{3}{26}\)
- 因為\(P(R,F)=P(Red)\times P(Face)=\frac{3}{26}\)，所以抽到紅色的牌與抽到有臉的牌是獨立事件。
1.8 邊際機率（Marginal Probability）
定義：在兩個或兩個以上的樣本空間中，只考慮某一個條件成立所發生的機率。
例子： 有一個樣本空間 \(S=\{1,2,3,4, 5,6,7,8,9,10\}\). \(P(A)=P(x=\rm {even})=0.5\). \(P(B)=P\{x\geq 7\}=0.4\)
使用councilor這筆資料，交叉分析年度與發包單位，計算邊際機率：
cr<-here::here('data','councilor.csv')
councilor<-read.csv(cr, sep=',',
header=TRUE, fileEncoding = 'BIG5')
tu<-table(councilor$Year, councilor$unit)
tu
##
## 公園處 水利處 新建工程處
## 2015 1 1 0
## 2016 0 0 8
margin.table(tu,1)/sum(tu)
##
## 2015 2016
## 0.2 0.8計算得知，2015年與2016年的邊際機率分別為0.2及0.8。
1.9 條件機率（Conditional Probability）
- 隨機變數Y出現y而另一個變數X出現x的機率，或者當事件B發生時事件A也發生，表示為\(P(A|B)\). \(A\) 發生是以\(B\)發生為前提。
- 條件機率可以用聯合機率以及邊際機率表示：
\[P(A|B)=\frac{P(A\cap B)}{P(B)}\]
- 從上一個式子可以推導出：
\[P(A\cap B)=P(A|B)\cdot P(B)=P(B\cap A)=P(B|A)\cdot P(A)\]
- 換句話說，條件機率可以是另一個條件機率乘以邊際機率再除以邊際機率，而這就是貝氏定理的基礎：
\[\begin{align} P(A|B)\cdot P(B) = P(B|A)\cdot P(A) \\ P(A|B)=\frac{P(B|A)\cdot P(A)}{P(B)} \end{align}\]
1.10 互賴事件 (Dependent Events)
例如當兵抽單位，前一個人抽到輕鬆單位，後面一個人抽到同樣單位的機率一定會變小，所以這兩個事件就是互賴事件，因為機率不相等。
例如台北市服務業集中，新北市有許多製造業，而男生又比較多從事製造業，女生比較多從事服務業，因此性別與職業是互賴事件，也就是男生這個事件的機率，會影響是不是從事製造業或者服務業的機率。
例如:從52張撲克牌中以抽出不放回的方式抽出兩張牌，請問第1張牌是皇后，跟第2張牌是皇后的事件是否為獨立事件？還是互賴事件？
假設\(P(E)\)代表第2張牌是皇后的機率，而\(P(A)\)是第一張牌為皇后的機率，因此\(P(A')\)是第一張牌不是皇后的機率。
已知\(P(A)=\frac{4}{52}\), 且\(P(A')=1-P(A)=\frac{48}{52}\), \(P(E|A)=\frac{3}{51}\), \(P(E|A')=\frac{4}{51}\)。
因此：
\[\begin{align} P(E) & = P(E\cap A)+P(E\cap A')\\ & = P(A)\cdot P(E|A)+P(B)\cdot P(E|B) \\ & = \frac{4}{52}\cdot\frac{3}{51}+\frac{48}{52}\cdot\frac{4}{51}\\ & = \frac{4}{52}\cdot\frac{3+48}{51} \\ & = \frac{4}{52} \end{align}\]
- 結論：如果52張撲克牌第1張抽出不放回，第2張牌抽到皇后的機率是 \(\frac{4}{52}\)，所以第1張牌是皇后，跟第2張牌是皇后的事件是獨立事件。
2 隨機變數
\[\rm{AB}=\frac{Hit}{Hit+Out}\]
\[\frac{4}{(4+6)}=0.4\]
☛課間練習：請問降雨機率是什麼？
☛課間練習：請問你可以用籃球為例想像投籃命中率的變數嗎？
tab<-cbind(結果=c("(T,T)","(H,T)","(T,H)","(H,H)"),
正面次數=c(0,1,1,2),
相對次數=c(0.25,0.25,0.25,0.25))
knitr::kable(tab,format = 'simple',
booktabs=T,caption="兩枚硬幣投擲結果與正面次數")| 結果 | 正面次數 | 相對次數 |
|---|---|---|
| (T,T) | 0 | 0.25 |
| (H,T) | 1 | 0.25 |
| (T,H) | 1 | 0.25 |
| (H,H) | 2 | 0.25 |
3 間斷變數
3.1 間斷變數的期望值與變異數
3.1.1 期望值
以表3.1為例，計算丟擲兩枚硬幣出現正面的期望值。計算結果應為1。也就是說，我們預期出現正面的次數為1。
| 正面次數 | 相對次數或機率函數 | 相乘 |
|---|---|---|
| 0 | 0.25 | 0 |
| 1 | 0.5 | 0.5 |
| 2 | 0.25 | 0.5 |
| sum | 1 | 1 |
R模擬從兩個集合\(\{0,1\}\)各抽樣1000次，然後加總抽樣結果，並且畫成圖 3.1 如下：
a<-c(0,1)
sample.1 <- 0; sample.2 <- 0
Sum <- 0
set.seed(02138)
for (i in 1:1000){
sample.1[i]=sample(a,1)
}
set.seed(02139)
for (i in 1:1000){
sample.2[i]=sample(a,1)
}
Sum <- sample.1 + sample.2
print(mean(Sum))
## [1] 1.031
table(Sum)
## Sum
## 0 1 2
## 232 505 263
hist(Sum)
Figure 3.1: 丟兩枚硬幣1000次的正面次數長條圖一
ggplot2畫圖一次。
D<-data.frame(Sum)
ggplot(data=D, aes(x=Sum))+
stat_count(geom='bar', fill='#ff22ee')+
geom_text(aes(label=(..prop..)), stat='count', colour='white', vjust=1.6) +
scale_y_continuous("",breaks=c(100,200,300,400,500),
labels = c(0.1,0.2,0.3,0.4,0.5))
Figure 3.2: 擲兩枚硬幣1000次的正面次數長條圖二
im1<-here::here('Fig','IMG_3167.JPG')
knitr::include_graphics(im1)
Figure 3.3: 打珠台
im2<-here::here('Fig','IMG_3168.JPG')
knitr::include_graphics(im2)
Figure 3.4: 打珠台事件與獎品
3.1.2 變異數
\[\begin{align} \boxed{V(X)=\sum_{i=1}^n (x_{i}-\mu)^2f(x_{i})} \end{align}\]
- 上面的變異數公式可以改寫為：
\[\begin{align} V(X) & =\sum_{i=1}^n (x_{i}-\mu)^2f(x_{i})\nonumber \\ & = \sum_{i=1}^n x_{i}^2f(x_{i})-\mu^2 \end{align}\]
或者
- 我們需要計算\(x_{i}^2f(x_{i})\)，就能計算變異數，以表3.2表示：
| 正面次數 | 相對次數或機率函數 | 平方項相乘 |
|---|---|---|
| 0 | 0.25 | 0 |
| 1 | 0.5 | 0.5 |
| 2 | 0.25 | 1 |
| sum | 1 | 1.5 |
- \(x_{i}^2f(x_{i})\)=1.5，而平均數等於1，所以變異數為1.5-1=0.5。
- 如果用變異數計算公式，則為：
\[\begin{align*} V(X) & =\sum_{i=1}^n (x_{i}-\mu)^2f(x_{i}) \\ & = \sum_{i=1}^n x_{i}^2f(x_{i})-\mu^2 \\ & = (0+0.5+1)-1 \\ & = 0.5 \end{align*}\]
- 我們可以用
R模擬的結果驗證，計算上述抽樣得到的結果為：
var(Sum)[1] 0.4945
\[f(x)=P(X=x)\] \[0\le P(X=x)\le 1\]
而且
\[\sum_{i}f(x_{i})=1\]
- 間斷變數常用的機率分配有：二元（白努利）分佈、二項分佈、Poisson分佈、超幾何分佈(hypergeometric)等等。
★ 二元（白努利）分佈(Bernoulli distribution)
- 白努利分佈假設成功事件的值X為1，其他為0。假設\(P(X=1)=p\)，那麼\(P(X=0)=1-p\)，PMF寫成：
\[\begin{align*} p(X=x)=p^x(1-p)^{1-x} = \begin{cases} p & \text{if}\quad x=1 \\ 1-p & \text{if}\quad x=0 \end{cases} \end{align*}\]
- \(x\)只有0, 1。當\(x=1\)，\(p(X=1)=p\)。當\(x=0\)，\(p(X=0)=1-p\)。我們可以求出二元分布的平均值為\(p\)，因為\(1\times p+0\times (1-p)=p\)。變異數為\(p(1-p)\)。
★ 二項分佈(Binomial distribution)
- 假如已知某個事件為真的機率為\(p\)，我們進行\(n\)次抽樣，其中\(X=k\)次得到某事件為真，該隨機變數服從二項分布。二項分布的機率質量函數寫成：
\[p(X=k)= \binom{n}{k}p^k(1-p)^{n-k}=\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}\]
- 假設酒駕致死的機率為\(p=82\%\)，如果警察隨機攔檢車輛發現有20位酒駕者，且假設攔檢結果彼此互相獨立，剛好出現12件致死車禍的機率是多少？
答：20次當中出現12次酒駕的機率質量函數為\[P(X=12)=\binom{20}{12}0.82^{12}(1-0.82)^{20-12}\]
可用R的
dbinom(12, 20, 0.82)
## [1] 0.01283- 可以畫圖表示如圖 3.5，次數為12與20的機率最低，17最高。
n=20; p=.82; x=12:n; prob=dbinom(x,n,p);
barplot(prob,names.arg = x,main="Binomial Barplot\n(n=20, p=0.82)",col="#ab3210")
Figure 3.5: 二項分配之次數與機率圖
- 也可以直接用指令計算，注意\(\frac{n!}{(n-k)!}\)的指令為
choose(n,k) ：
p=0.82
n=20
k=12
choose(n, k)*(p^k)*(1-p)^(n-k)[1] 0.01283
#or
factorial(n)/(factorial(k)*factorial(n-k))*(p^k)*(1-p)^(n-k)[1] 0.01283
- 二項分布的期望值剛好是二元分佈的期望值的\(n\)倍，\(n\)是實驗的次數，所以寫成\(E(X)=np\)。而變異數\(V(X)=np(1-p)\)。
★ Poisson分佈
Poisson是一段期間內發生事件的次數的機率分佈，事件非同時發生，且發生或者不發生是二元的變數，所以可視作白努利試驗，n個事件的組合，也可以稱做二項式分布。例如從9點到10點到郵局的人數，晚上11點到12點登入PTT的人數。
Poisson分佈的PMF寫成：
\[p(x=k;\lambda)=\frac{e^{-\lambda}\lambda^{k}}{k!}\quad \rm{for}\quad k=0,1,2,\ldots\]
- 其中Poisson分佈的平均值為\(\lambda\)，\(\lambda\)是在一定期間內的事件發生次數的平均數。圖3.6顯示三種平均數的Poisson分佈：
N=1000
d=tibble(lambda=seq(1,18, by=6))
random_poisson=function(lambda,n) {
rpois(n,lambda)
}
d = d %>% mutate(randoms=purrr::map(lambda,random_poisson, N))
#d %>% purrr::unnest(randoms)
d %>% tidyr::unnest(randoms) %>%
ggplot(aes(x=randoms))+
geom_histogram(aes(fill=as.factor(lambda)))+
facet_wrap(~lambda,scales="free") +
theme(legend.position = 'none')
Figure 3.6: 3種平均數的Poisson分佈
平均數\(\lambda\)越大，越接近常態分佈。
現代統計學 第197頁，已知平均5分鐘有1.5人使用ATM，半小時平均有9人，請問半小時有5人使用ATM的機率是多少？我們可以應用dpois(k, \(\lambda\)) ，求出當平均值為9的時候，x=5的機率：
lambda=9
dpois(5, lambda)[1] 0.06073
- 或者用公式計算：
lambda=9
k=5
p=lambda^{k}*exp(-lambda)/factorial(k)
p[1] 0.06073
現代統計學 第202頁，已知咖啡廳平均每小時有50位顧客，請問每小時需要幾位服務人員，才不會讓每位服務人員需要服務超過8位顧客的機率高於0.2？換句話說，在\(k>8\)的條件下，\(\lambda\)應該等於多少，才不會出現一個服務生必須服務超過8個客人？而且\(\lambda=50/b\)，\(b\)是服務生人數。圖 3.7 顯示Poisson的機率分佈，也就是\(P(X=k, \lambda)\)。把k對應的機率累積，就可以得到圖 3.8 顯示的Poisson分佈的累積機率：\(f(X)=P(X>k|\lambda)\)。當\(k\)越大，\(P(X>k|\lambda)\)也越大，然後就接近1。曲線上的每一個點代表當X=k時，函數所累積的機率，例如當\(P(X>10|\lambda=10)=0.583\)，也就是\(P(X=0,\lambda=10)+P(X=1,\lambda=10)+\ldots+P(X=10,\lambda=10)\)。
k<-seq(0, 20, by=1)
y <- dpois(k, lambda=10)
plot(k, y, type='l', lwd=1.5, col='#1a00ab',
ylab=expression(paste('P(X=k|', lambda,')')))
Figure 3.7: Poisson機率分佈之機率密度圖
k<-seq(0, 20, by=1)
y <- ppois(k, lambda=10)
plot(k, y, type='l', lwd=1.5, col='#1a00ab',
ylab=expression(paste('P(X>k|', lambda,')')))
abline(v=10, lty=2)
abline(h=0.583, lty=2)
Figure 3.8: Poisson機率分佈之累積機率密度圖
- 回到問題，我們可以應用
ppois(k=8, \(\lambda\)) ，求出當平均值為\(50/y\)、k>8的時候，y應該等於多少，才會讓顧客被服務的平均人數至少為8的機率不到0.2。這裡要注意的是，ppois(k=8, \(\lambda\)) 是從0累加到8，也就是顧客被服務的平均人數0一直到8，平均人數越多應該需要越多服務人員。但其實我們關心的是超過8之後的機率是否小於0.2，因此要用1來減去累積機率，或者直接要R計算從曲線另一邊累積的機率。圖3.9顯示當b=8，\(1-P(8;6.25)\)的累積機率小於0.2。當然，咖啡廳可以雇用更多服務人員，\(1-P(8;6.25)\)的累積機率也越小，也就是不太可能出現一位服務人員需要招呼超過8為客人的情況，但是可能會超出雇主的負擔。
b <- seq (0, 10, by = 1)
y<-c()
for (i in 1:11){
y[i] <- ppois(b[i], lambda=c(50/b[i]), lower.tail=F)
}
plot(b, y, type='l', lwd=1.5, col='#1a00ab',
ylab=expression(paste('1-P(X>k|', lambda,')')))
abline(h=0.2, lty=2)
abline(v=8, lty=2)
Figure 3.9: Poisson機率分佈之累積機率密度圖
- 當\(n\geq 20\)時，Poisson分佈可以代替二項分佈。例如當機率等於0.45，擲100次硬幣得到10次正面的機率為0.048，而以Poisson機率分佈而言，當\(\lambda=45\)，事件發生40次的機率則是0.047。這是因為二項分佈的期待值為\(np=100*0.45=45\)，所以當平均值接近、n越大（隨機實驗次數越多），兩個機率分佈越接近。
dbinom(40, size =100, p=0.45)
## [1] 0.0488
dpois(40, 45)
## [1] 0.047164 連續變數
身高、體重、體溫、收入、雨量、時間、距離等等是常見的連續變數，因為這些變數的個數是無限而且不可數。
連續變數常用的機率分配有：常態分佈、標準常態分佈、均等分佈、指數分佈等等。
定義：連續隨機變數\(X\)，其值落在\([a, b]\)，也就是\(a\leq X\leq b\)。如果此函數非負而且在區間\([a, b]\)連續，而且如果 \(\int_a^b f(x)dx=1\)，則稱\(f(x)\)為\(X\)的機率密度函數(probability density function, pdf)。
用\(X\)代表變數，\(X\)可能是某一個整數(\(0,1,2,\ldots,100\))，或者無理數(\(\frac{1}{\sqrt{10}}\))，我們想像每一個變數的值有若干機率，假設有\(101\)個值，機率可以寫成
\[\rm{Pr}(X=i)=\frac{1}{101}\quad for \hspace{0.4em}i=0,1,\ldots,100\]
- 如果\(x\)超過100，那麼變數\(X\)不可能是\(x\)，換句話說，
\[\rm{Pr}(X=x)=0\quad if\hspace{.4em} x\hspace{.4em} is\hspace{.4em} not\hspace{.4em} in\hspace{.4em} \{1,\ldots,n\}\]
如果\(x\in(0,1)\)，也就是在0跟1之間，也就是說\(0\le x\le 1\)，那麼當\(x>1\)，\(Pr(X=x)=0\)；當\(x<0\)，\(Pr(X=x)=0\)。
不僅如此，我們可以說\(Pr(X=x)=0\)，這是因為在連續變數中，我們很難找到\(x\)剛好等於某一個實數，或者說一條直線的面積是0。但是我們可以用積分計算包含一個區間的面積，也就是\(Pr(a\le X\le b)=\int_a^bf(x)dx\)。
\(f(x)\)是一個機率分佈(Probability distribution)，而不是固定的數字，\(f(x)\)與\(X\)這兩者之間的關係稱為機率密度(probability density)。進一步假設\(x\)位於一個區間\(I\)，這個區間的大小以及密度\(f(x)\)決定變數\(X\)的值落在這個區間的機率，表示成：
\[\rm{Pr}(X\in I)\approx f(x)\times \rm{Length\hspace{.5em} of \hspace{.5em}I}\]
- 例如常態（高斯）分佈的機率密度函數是：
\[f(x)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x-\mu)^2}{2\sigma^2})\]
- 隨機變數如果是常態分佈，寫成：
\[X\sim N(\mu, \sigma^2)\]
- 我們可以透過\(W=\frac{x-\mu}{\sigma}\)轉換常態分佈為標準常態分佈，\(X\sim N(0, 1)\)也就是平均值為0、變異數為1的常態分佈(轉換過程)，寫成：
\[\phi(z)=\frac{1}{\sqrt{2\pi}}exp(-\frac{1}{2}z^2)\]
- 我們可以寫
R表示標準常態分佈的機率分佈密度函數為：
s=1; mu=0
f <- function(x){
1/sqrt(2*pi*s^2)*exp(-(1/2)*(x-mu)^2/s^2)
}- 然後用
integrate() 計算積分，可得1，表示在負無限大與無限大之間，它是z的機率分佈密度函數。
integrate(f, lower=-Inf, upper=Inf)
## 1 with absolute error < 9.4e-05- 注意並不是只有標準常態分佈，曲線下的面積才會等於1，例如\(X\sim N(0, 0.5)\)，積分PDF結果為1：
integrate(dnorm, 0, 1, lower=-Inf, upper=Inf)
## 1 with absolute error < 9.4e-05
integrate(dnorm, 0, 0.5, lower=-Inf, upper=Inf)
## 1 with absolute error < 3.1e-05
integrate(dnorm, 2, 0.5, lower=-Inf, upper=Inf)
## 1 with absolute error < 2e-06- 可以換其他的變異數，得到一樣的答案：
mu<-0; s=0.5
f <- function(x){
1/sqrt(2*pi*s^2)*exp(-(1/2)*(x-mu)^2/s^2)
}
integrate(f, lower=-Inf, upper=Inf)
## 1 with absolute error < 3.1e-05
mu<-2; s=0.5
integrate(f, lower=-Inf, upper=Inf)
## 1 with absolute error < 2e-06dnorm() 產生常態分佈的機率密度函數\(f(x)\)，我們可以用在積分：
integrate(dnorm, 0, 1, lower=-Inf, upper=Inf)
## 1 with absolute error < 9.4e-05
integrate(dnorm, 0, 0.5, lower=-Inf, upper=Inf)
## 1 with absolute error < 3.1e-05
integrate(dnorm, 2, 0.5, lower=-Inf, upper=Inf)
## 1 with absolute error < 2e-06☛如果\(x\in(1,\infty)\)，請問以下這個函數是連續函數嗎？
\[h(x)=\frac{3}{x^4}\]
- 我們從常態分佈中抽出1,000個值，然後畫長條圖如圖 4.1。注意兩個圖有不同的縱軸，圖一是相對次數，圖二是相對次數密度，相對次數=相對次數密度乘以組距：
set.seed(2020)
extmp<-rnorm(mean=10, sd=1, 1000)
par(mfrow=c(1,2))
barplot(table(cut(extmp, 15))/1000, ylab="相對次數",main="圖1", family="YouYuan")
hist(extmp, main="圖2", freq = F, ylab="相對次數密度", family="YouYuan")
Figure 4.1: 單位與實數組距的長條圖
- 假設我們相信\(X\)可能落在某一個範圍內，例如稱作\(A\)，\(A\)從a, b，也就是把a到b的區域的\(f(x)\)加總起來：
\[\rm{Pr}(x \in A) =\int_{A}p(x)dx=\int_{a}^{b}p(x)dx\]
- 例如常態分佈的部分面積如圖 4.2：
par(family="YouYuan" )
#x <- seq(-3, 3,0.1)
curve( dnorm(x),
xlim = c(-3, 3),
ylab = "Density",
main = "機率密度與區域",
col='red', lwd=2)
cord.1x <- c(0,seq(0,1,0.01),1)
cord.1y <- c(0,dnorm(seq(0,1,0.01)),0)
# Make a curve
#curve(dnorm(x,0,1), xlim=c(-3,3), main='Standard Normal',lwd=2)
# Add the shaded area.
polygon(cord.1x,cord.1y,col='#44aa00ee')
Figure 4.2: 機率分佈曲線下的區域
- 我們可以用積分求0到1之間曲線底下的面積，得到\(0\le x\le 1\)的機率，也就是套用積分在常態分佈的函數：
\[f(x)=\frac{1}{\sqrt{2\pi}\sigma}exp(-\frac{(x-\mu)^2}{2\sigma^2})\]
mu<-0; s=1
f <- function(x){
1/sqrt(2*pi*s^2)*exp(-(1/2)*(x-mu)^2/s^2)
}
integrate(f, lower=0, upper=1)
## 0.3413 with absolute error < 3.8e-15結果約為0.34，也就是說在0與1個標準差之間有\(34\%\)的面積，如果是0左右1個標準差，那麼面積就是\(68\%\)。
累加機率函數(cumulative density function, CDF)代表機率從F(a)\(=0\)一直累加，最後到F(b)\(=1\)。表示為：
\[F(X)=P(X\leq x)=\int_a^xf(t)dt\]
如果c<d，則F(c\(\le\)X\(\le\)d)=F(d)-F(c)
F(x)=P(X\(\le\)x)
對cdf進行微分可得機率密度函數表示為：
\[f(x)=\frac{d}{dx}\int_{-\infty}^xf(x)dx\]
- 期望值：
\[E(X)=\int_a^bxf(x)dx=\mu\]
- 變異數：
\[\begin{align} V(X) & =\int_a^b(x-\mu)\cdot f(x)dx \\ &=\int_a^bx^2\cdot f(x)dx-\mu^2 \\ &=E(X^2)-[E(X)]^2 \end{align}\]
- 參考林惠玲與陳正倉(2004)，統計學，假設公車每10分鐘來一班，也就是每一分鐘有\(\frac{1}{10}\)的機會等到公車，機率密度函數表示為：
\[\begin{equation*} f(x)= \begin{cases} \frac{1}{10} & 0\le x\le 10\\ 0 & \text{otherwise} \end{cases} \end{equation*}\]
- 等公車的時間不超過3分鐘的機率可視為累加不超過3分鐘的機率：
\[\begin{align} F(X)=P(X\le 3)& =\int_a^xf(x)dx \\ & = \int_a^x \frac{1}{10}dx \\ & = \frac{x}{10} \\ & = \frac{3}{10} \end{align}\]
- 等公車的時間超過4分鐘的機率可視為1扣掉累加不超過4分鐘的機率：
\[\begin{align} 1-F(X=4) & = 1-P(X\le 4) \\ & = 1 - \frac{4}{10} \\ & = 0.6 \end{align}\]
- 平均等待公車的時間為多少？
\[\begin{align} E(X) & = \int_a^bxf(x)dx \\ & = \int_{0}^{10}x\frac{1}{10}dx \\ & = \frac{1}{2}\cdot \frac{1}{10}x^2 \\ & = \frac{100}{20}\\ & = 5 \end{align}\]
- 等待公車的時間變異數為多少？
\[\begin{align} V(X) & = \int_a^b (x-\mu)^2\cdot f(x)dx \\ & = \int_{0}^{10} (x-5)^2\cdot \frac{1}{10}dx \\ & = \int_{0}^{10} (x^2-10x+25)\cdot \frac{1}{10}dx \\ & = \frac{1}{30}x^3-\frac{1}{2}x^2+2.5x \\ & = \frac{1000}{30} -\frac{100}{2}+25 \\ & = 8.33 \end{align}\]
- 當我們對\(f(x)=a+bx+cx^2\)微分，寫成\(f(x')\)或者\(\frac{d}{dx}f(x)\)，方法為把x降1次方，也就是\(x^1\)降為\(x^{1-1}=x^0=1\)。同時，常數視為0：
\[\begin{align} f(x') & = 0+ bx^{1-1} + cx^{2-1} \\ & = b+2cx \end{align}\]
- 反過來，當我們對\(f(x)=a+bx+cx^2\)積分，方法為把x加1次方，也就是\(x^1\)變成為\(x^{1+1}=x^2\)，同時，x前面要乘以\(\frac{1}{1+1}=\frac{1}{2}\)。常數則是幫它加上x的1次方：
\[\begin{align} \int_{-\infty}^{\infty} (a+bx+cx^2)dx & = ax+\frac{1}{2}bx^2+\frac{1}{3}cx^3 \end{align}\]
4.1 連續變數的期望值與變異數
4.1.1 期望值
\[\boxed{E(X)=\int_{a}^{b}xf(x)dx=\mu} \]
- 在間斷變數的例子，變數有個別的變量，例如變數\(X\)有四個可能的值1, 2, 3, 4，發生的機率如表 4.1：
| X | 1.0 | 2.0 | 3.0 | 4.0 |
| p | 0.1 | 0.3 | 0.4 | 0.2 |
tmp <- c(rep(1,1), rep(2, 3), rep(3,4), rep(4,2))
ttmp <- table(tmp/10)
names(ttmp)<-c(1,2,3,4)
barplot(ttmp, col='#cc002233')
Figure 4.3: 間斷變數之相對次數圖
\[\begin{equation*} \begin{split} p(2\leq t\leq 4) & = \frac{1}{36}\int_2^4(-x^2+6x)\, dx \\ & = \frac{1}{36}(-\frac{4^3}{3}+\frac{1}{2}(6\cdot 4^2)-(-\frac{2^3}{3}+\frac{1}{2}(6\cdot 2^2)))\\ & = \frac{1}{36}(-\frac{64}{3}+48+\frac{8}{3}-12)\\ &=\frac{13}{27} \\ &\approx 0.481 \end{split} \end{equation*}\]
- 圖 4.4 顯示四種大小不同的\(x\)變量：
par(mfrow=c(2,2))
x<-seq(from=-2, to=6, length.out = 1000)
y=(1/36)*(-(x^2)+6*x)
plot(x, y, ylim=c(0, 0.3), lwd=1, cex=0.1)
x<-seq(from=0, to=6, length.out = 1000)
y=(1/36)*(-(x^2)+6*x)
plot(x, y, ylim=c(0, 0.3), lwd=1, cex=0.1)
x<-seq(from=1, to=6, length.out = 1000)
y=(1/36)*(-(x^2)+6*x)
plot(x, y, ylim=c(0, 0.3), lwd=1, cex=0.1)
x<-seq(from=3, to=6, length.out = 1000)
y=(1/36)*(-(x^2)+6*x)
plot(x, y, ylim=c(0, 0.3), lwd=1, cex=0.1)
Figure 4.4: 連續變數之機率密度圖
x<-seq(from=-1, to=6, length.out = 1000)
y=(1/36)*(-(x^2)+6*x)
plot(x, y, ylim=c(0, 0.3), lwd=1, cex=0.1)
cord.1x <- c(2, seq(2,4,0.01), 4)
cord.1y <- (1/36)*(-(cord.1x^2)+6*cord.1x)
cord.1y[1]<-0; cord.1y[203]<-0
polygon(cord.1x,cord.1y,col='#aa11cc33')
Figure 4.5: 機率密度圖的區域
R提供積分以及微分的功能，在此我們用積分計算\(\int_{2}^{4}f(x)dx\)，也就是\(F(4)-F(2)\)如下：
F <- function(x) (1/36)*(-(x^2)+6*x)
integrate(F, 2, 4)
## 0.4815 with absolute error < 5.3e-15- 得到的結果跟我們自行計算一樣為0.481。也就是說，\(2\le x\le4\)的機率為\(48.1\%\)，接近五成。要注意的是，期望值公式為
\[\boxed{E(X)=\int_{a}^{b}xf(x)dx=\mu} \]
R計算。
4.1.2 變異數
- 在了解連續變數的平均數之後，我們想要知道變數的離散程度。連續變數的變異數的計算方式為：
\[\begin{align} V(X)=\int_{a}^{b}(x-\mu)^2\cdot f(x)dx \tag{4.1} \end{align}\]
- 方程式 (4.1)可以改為：
\[\begin{align*} V(X) & =\int_{a}^{b}(x-\mu)^2\cdot f(x)dx \\ & = \int_{a}^{b}x^2\cdot f(x)dx-\mu^2 \\ & = E(X^2) -[E(X)]^2\tag{4.2} \end{align*}\]
R內建的變異數公式對照：
set.seed(02138)
x<-rnorm(100, 0, 2)
xbar<-mean(x)
Meansquare<-sum(x^2)/100
squareofmean<-xbar^2
cat("Variance=",Meansquare - squareofmean,"\n")
## Variance= 3.948
cat("Variance of sample=",var(x))
## Variance of sample= 3.988R的樣本變異數公式為：\(\frac{\sum (X-\bar{X})^2}{n-1}\)，所以跟母體變異數公式\(\frac{\sum (X-\bar{X})^2}{n}\)得到的結果略有差異。
- 注意\(E(X^2)\)指的是累加\(X_1^2,X_2^2,\ldots,X_{i}^2\)然後除以\(n\)，所以：
\[E(X^2)=\int x^2f(x)dx\]
\(\blacksquare\)假設有一個連續函數如下：
\[\begin{equation*} f(x) = \begin{cases} 2(1-x) & \text{if}\hspace{1em} 0\le x\le 1\\ 0 & \text{otherwise} \end{cases} \end{equation*}\]
- 首先我們先確認\(\int_{0}^{1}f(x)dx=1\)
f<-function(x)2*(1-x)
integrate(f, 0, 1)1 with absolute error < 1.1e-14
- 然後根據期望值公式：
\[\boxed{E(X)=\int_{a}^{b}xf(x)dx=\mu} \]
\[\begin{align*} E(X) & =\int_{0}^{1}x[2(1-x)]dx \\ & = 2\int_{0}^{1}(x-x^2)dx \\ & = 2[(\frac{1^2}{2}-\frac{1^3}{3})- (\frac{0^2}{2}-\frac{0^3}{3})] \\ & = 1/3 \end{align*}\]
- 接下來計算變異數，公式為：
\[\boxed{V(X) =\int_{-\infty}^{\infty}(x-E(X))^2\cdot f(x)dx}\]
- 因為平均數等於\(\frac{1}{3}\)，我們代入\(E(X)=\frac{1}{3}\)，而\(f(x)=2(1-x)\)，求變異數如下：
\[\begin{align*} V(X) & = \int_{0}^{1}(x-\frac{1}{3})^2\cdot 2(1-x)dx \\ & = 2\int_{0}^{1}(x^2-\frac{2}{3}x+\frac{1}{9})\cdot (1-x)dx \\ & = 2(-\frac{1}{4}x^4+\frac{5}{9}x^3-\frac{7}{18}x^2+\frac{1}{9}x)\Biggr|_{0}^{1} \\ & = \frac{1}{18} \end{align*}\]
- 用
integrate 函式驗證，得到0.055，約為\(\frac{1}{18}\)。
f<-function(x)((x-1/3)^2)*2*(1-x)
integrate(f, lower=0, upper=1)0.05556 with absolute error < 6.2e-16
☛假設有一個連續函數如下：
\[\begin{equation*} f(Y) = \begin{cases} \frac{3}{Y^4} & \text{if}\hspace{1em} x\ge 1\\ 0 & \text{otherwise} \end{cases} \end{equation*}\]
- 請求出\(f(Y)\)的期望值以及變異數？
4.2 常態機率分佈
- 常態分佈的值落在平均數加減3個標準差等距的範圍有0.9974的機率，表示為：
\[P(\mu-3\sigma<X<\mu+3\sigma)=0.9974\]
- 常態分佈的值落在平均數加減2個標準差等距的範圍有0.9544的機率，表示為：
\[P(\mu-2\sigma<X<\mu+2\sigma)=0.9544\]
- 常態分佈的值落在平均數加減1個標準差等距的範圍有0.6826的機率，表示為：
\[P(\mu-\sigma<X<\mu+\sigma)=0.6826\]
x<-seq(from=-3, to=3,length.out=1000)
s=1
y <- (1/sqrt(2*pi)*s)*(exp(-x^2/2*s^2))
ylim<-c(0, 0.85)
plot(x, y ,main=" ", cex=0.1, ylim=ylim)
s=1.4
curve((1/sqrt(2*pi)*s)*(exp(-x^2/2*s^2)), from=-4, to=4, col="#CC2233FF", add=T)
s=2
curve((1/sqrt(2*pi)*s)*(exp(-x^2/2*s^2)), from=-3, to=3, col="blue", add=T)
Figure 4.6: 常態機率分佈
- 如果資料近似鐘形的常態分佈，可以在長條圖上面加上常態分佈曲線，如圖 4.7：
#抽樣然後加總得到2000個數字
set.seed(02138)
ds <-c()
for (i in 1:2000){
ds[i] <- sum(sample(1:6, 1), sample(1:6,1),
sample(1:6,1), sample(1:6, 1))
}
#轉成資料框
gnorm<-data.frame(dice=ds)
#畫長條圖，注意y是機率密度
ggplot2::ggplot(aes(x=dice, y=..density..), data=gnorm) +
geom_histogram(bins=20, fill='lightblue') +
geom_density(aes(y=..density..), col='steelblue') +
theme_bw()
Figure 4.7: 近似常態分佈之長條圖加上常態分佈曲線
x <-seq(from=-3, to=3, 0.01)
y<-c()
f.normal <- function(x) {
y <- 1/(sqrt(2 * pi)) * exp(-0.5 * x^2)
}
plot(x, cumsum(f.normal(x))/100, cex=0.2)
Figure 4.8: 累積機率密度圖
R提供的函數畫累積密度曲線，如圖 4.9。當\(x=0\)，也就是平均數，累積密度等於0.5。
curve(pnorm, from=-3, to=3, col='#EE11CCFF', lwd=1.5)
segments(0, -0.03, 0, 0.5, lty=2, lwd=1.5)
segments(-3.5, 0.5, 0, 0.5, lty=2, lwd=1.5)
Figure 4.9: R函數之累積密度曲線圖
- 圖 4.10 顯示當x=0時，面積剛好等於0.5。
#x <- seq(-3, 3,0.1)
curve( dnorm(x),
xlim = c(-3, 3),
ylab = "Density",
#main = "機率密度與區域",
col='red', lwd=2)
cord.1x <- c(0,seq(-3,0,0.01),0)
cord.1y <- c(0,dnorm(seq(-3,0,0.01)),0)
# Add the shaded area.
polygon(cord.1x,cord.1y,col='#33aacc')
Figure 4.10: 機率分佈曲線下的一半區域
4.2.1 常態分配相關函式
R的常態分佈函式有dnorm 產生機率密度函數，也就是代入X之後的值。我們以IQ這個連續變數為例，最小為70、最高為150，平均值110、標準差則是任意的15。我們可以產生對應的機率密度，並且畫成圖4.11：
sample.range <- seq(70, 150, by=1)
iq.mean <- 110
iq.sd <- 15
iq.dist <- dnorm(sample.range, mean = iq.mean, sd = iq.sd)
iq.df <- data.frame("IQ" = sample.range, "Density" = iq.dist)
ggplot2::ggplot(iq.df, aes(x = IQ, y = Density)) + geom_line(colour="#21ab01") 
Figure 4.11: IQ機率密度圖
pnorm 產生累積機率密度函數，也就是代入X之後產生從最小值一直累積到X的機率。畫成圖 4.12：
cdf <- pnorm(sample.range, iq.mean, iq.sd)
iq.df <- cbind(iq.df, "CDF_LowerTail" = cdf)
ggplot(iq.df, aes(x = IQ, y = CDF_LowerTail)) + geom_line(colour="#21ab01") 
Figure 4.12: IQ累積機率密度圖
pnorm 產生累積機率密度函數，qnorm 則是pnorm 的相反，也就是產生從最小值一直累積到X的機率所對應的值。畫成圖4.13：
prob.range <- seq(0, 1, 0.001)
icdf.df <- data.frame("Probability" = prob.range, "IQ" = qnorm(prob.range, iq.mean, iq.sd))
ggplot(icdf.df, aes(x = Probability, y = IQ)) + geom_line(color="#bb012e")
Figure 4.13: IQ累積機率密度圖
- 為了確認
pnorm 與qnorm 的結果，我們計算第25百分位的IQ，然後反過來計算該IQ值對應的累積百分比如下：
# the 25th IQ percentile?
print(icdf.df$IQ[icdf.df$Probability == 0.25])
## [1] 99.88
# x=? 25th IQ percentile?
perc <- function (x){
paste0(round(x * 100, 3), "%")
}
perc(iq.df$CDF_LowerTail[iq.df$IQ == 100])
## [1] "25.249%"可見得智商100大概落在25百分位，只贏25%的全部人。
rnorm 產生來自常態分佈的隨機亂數，我們產生三個不同樣本規模但是相同的平均數以及變異數的樣本，然後畫相對次數的長條圖 4.14，可以發現，樣本規模越大，則越集中在平均數附近。
set.seed(1)
n.samples <- c(100, 1000, 10000)
my.df <- do.call(rbind, lapply(n.samples, function(x) data.frame("SampleSize" = x, "IQ" = rnorm(x, iq.mean, iq.sd))))
# show one facet per random sample of a given size
ggplot(data = my.df, aes(x =IQ)) +
geom_histogram(aes(y =stat(count) / sum(count), fill=as.factor(SampleSize))) +
facet_wrap(.~SampleSize, scales = "free_y")+
scale_fill_discrete(name='')
Figure 4.14: 三種樣本規模的長條圖
4.3 均等機率分佈
- 均等分配指隨機變數在某一連續的區間，有同等的機率密度，例如等車時間、飛行時間、排隊買美食時間。可以寫成:
\[X\sim U(a, b)\]
- 例：某個月的日出時間是在6點半到7點之間，所以在這個區間的日出機率為\(\frac{1}{30}\)。如果6點50分才起床，看到日出的機率是\(P(X\leq 10)=\frac{10}{30}=\frac{1}{3}\)。
\[\begin{equation*} f(x) = \begin{cases} \frac{1}{b-a} & a\le x\le b\\ 0 & \text{otherwise} \end{cases}\tag{4.3} \end{equation*}\]
- 在方程式(4.3)中，如果a=0、b=1，\(f(x)=1\)。換句話說，在\(\{0,1\}\)這個區間內的\(X\)，對應到\(f(x)\)的值都等於 1。我們可用
dunif() 表示機率密度：
x<-seq(0,1, by=0.1)
dunif(x, min=0, max=1)
## [1] 1 1 1 1 1 1 1 1 1 1 1既然\(f(x)=1\)，那麼\(\int_{0}^{1}f(x)dx=1\times 1=1\)。證明了\(f(x)=\frac{1}{b-a}\)是一個連續隨機變數。
假設\(X\sim U(0,2)，\)圖 4.15 顯示均等機率分佈曲線下的區域：
x<-seq(from=-1,to=3,length.out=1000)
ylim<-c(0,0.6)
plot(x,dunif(x,min=0,max=2),main=" ",
type="l",ylim=ylim, yaxt='n',
ylab=expression(paste('f',(x))))
axis(side=2, at=seq(0.1,0.5,by=0.2),
labels = seq(0.1,0.5,by=0.2))
Figure 4.15: 均等機率分佈曲線下的區域
- 由以下可知，\(f(x)=0.5\quad \forall 0\leq X\leq 2\)
x=seq(0,2, by=0.1)
dunif(x,min=0,max=2)[1] 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 [20] 0.5 0.5
- 圖 4.16 顯示日出時間是6點半到7點之間，而在6點50分到7點之間的均等機率分佈曲線下的區域以綠色代表：
x<-seq(from=-2,to=32,length.out=1000)
ylim<-c(0,0.045)
plot(x,dunif(x,min=0,max=30),main=" ",type="l",
ylim=ylim)
cord.1x <- c(20,seq(20,30,0.01),30)
cord.1y <- c(0,dunif(seq(20,30,0.01),
min=0, max=30),0)
polygon(cord.1x,cord.1y,col='#ace100')
cord.1x <- c(0,seq(0,20,0.01),20)
cord.1y <- c(0,dunif(seq(0,20,0.01),
min=0, max=30),0)
polygon(cord.1x,cord.1y,col='#ee22cc00')
Figure 4.16: 一致機率分佈曲線下的區域
4.3.1 均等分配的期望值與變異數
- 均等分配的期望值為：
\[\boxed{E(X)=\frac{a+b}{2}}\]
\[\boxed{V(X)=\frac{(b-a)^2}{12}}\]
- 也可以寫成：
\[E(X^2)-E(X)^2\]
- 其中需要求出\(X^2\)的期望值如下：
\[\begin{align*} E(X^2)&=\int_{a}^{b}X^2f(x)dx \\ & = \int_{a}^{b}\frac{X^2}{b-a}dx \\ & = \frac{1}{b-a}\times\frac{1}{3}\int_{a}^{b}X^3 \\ & = \frac{1}{b-a}\times\frac{1}{3}(a-b)^3 \tag{4.4} \end{align*}\]
\[\begin{equation*} f(x) = \begin{cases} \frac{1}{X} & 0\le x\le 1\\ 0 & \text{otherwise} \end{cases}\tag{4.3} \end{equation*}\]
- 求出期待值？
因為\[\int_{0}^{1}zd(z)\] 所以 \[E(z)=\frac{0+1}{2}=0.5\]
- 求出變異數？
\[E(z)=\frac{(1-0)^2}{12}=0.083\]
- 或者根據方程式(4.4)：
\[\begin{align*} E(Z^2)-E(z)^2 & =\frac{1}{1-0}\frac{1}{3}(1-0)^2-0.5^2\\ & = 0.333-0.025 \\ & = 0.083 \end{align*}\]
g<-function(x){x}#f(x)
h<-function(x){x^2}#f(x^2)
eg<-integrate(g, 0, 1)
eh<-integrate(h, 0, 1)
eg; eh0.5 with absolute error < 5.6e-15 0.3333 with absolute error < 3.7e-15
#variance
as.numeric(eh[1])-as.numeric(eg[1])^2[1] 0.08333
4.3.2 均等分配相關函式
dunif(x, min, max) 產生機率密度函數，也就是代入X之後的值。例如公車一小時一班，在一個小時內任何一分鐘公車來的機率是：
dunif(1, min=0, max=60)[1] 0.01667
punif 回傳對應輸入的均勻分佈累積機率值，也就是\(P(X\leq x)\)，例如\(X\sim U(0,10)\)，則\(P(X\leq 1)=0.1\)，而\(P(X\leq 2)=0.2\)，以此類推：
px <- c(1, 2, 3)
fx <- punif(px, 0, 10)
df <- data.frame(px, fx)
ggplot2::ggplot(df, aes(x=px,y=fx)) +
geom_step() +
geom_point(data=df, mapping=aes(x=px, y=fx),
color="red") 
Figure 4.17: 均勻分佈的累積機率圖
qunif 回傳均勻分佈的百分位如圖4.18：
prob.range <- seq(0, 1, 0.001)
icdf.df <- data.frame("Probability" = prob.range, y = qunif(prob.range))
ggplot(icdf.df, aes(x = Probability, y = y)) + geom_line(color="#1dd02f")
Figure 4.18: 均勻分佈的百分位圖
runif 產生均等分配的隨機變數，例如抽出位於(0,10)的1000個數字，繪成長條圖4.19之後，每一個數字的機率密度約為0.1：
set.seed(1000)
rand.unif<-runif(1000, 0, 10)
hist(rand.unif, freq = FALSE, col='#4433ae', density = 40, xlab='x', main='')
Figure 4.19: 均等分配隨機變數的長條圖
punif(1, 0, 10)[1] 0.1
4.4 指數機率分佈
\[f(x)=\lambda e^{-\lambda x}\quad x\ge 0,\lambda>0\]
- 則\(f(x)\)為指數分配，其中\(\lambda\)為單位時間事件發生的平均數，(次數/單位時間)，也稱為速率參數(rate parameter)。但是\(\frac{1}{\lambda}\)是發生事件的平均時間。例如排班的計程車司機平均15分鐘可以載到一個乘客，\(\lambda=15\)。而對於司機來說，每分鐘平均可以載到\(\frac{1}{15}\)的乘客，\(\frac{1}{\lambda}=\frac{1}{15}\)。機率密度函數可以寫成：
\[f(x)=\frac{1}{15}e^{-\frac{x}{15}}\]
- 累積機率密度函數則是：
\[P(X\leq x_{0})=1-P(X> x_{0})=1-e^{-\lambda x_{0}}\]
- 指數分配的隨機變數寫成：
\(X\sim \text{Exp}(\lambda)\quad 0\leq X\leq \infty\)
4.4.1 指數分配期望值與變異數
\[E(X)=\frac{1}{\lambda}\]
\[V(X)=\frac{1}{\lambda^2}\]
4.4.2 指數分配的相關指令
dexp 回傳機率密度，以上一題為例，可假設\(0\leq X\leq 50\)，得到對應的機率密度由高而低如圖4.20：
x = seq(0, 50, by=0.1)
df<-data.table::data.table(x, y = dexp(x, rate=1/15))
ggplot(df, aes(x=x, y=y)) +
geom_line(color="#4013ae")
Figure 4.20: 指數分配機率密度圖
- 30分鐘內有客人上車的機率是累積機率\(P(X\leq x_{0}=30)=1-P(X>30)\)：
x=30; rate=1/15
pexp(x, rate)[1] 0.8647
- 公式計算：
rate=1/15; x=30
1-exp(-rate*x)[1] 0.8647
pexp 回傳累積機率密度，也就是\(F(X=x_{0})=P(X\leq x_{0})=1-P(X> x_{0})\)。以上一題為例，可假設\(0\leq X\leq 20\)，得到對應的機率密度由高而低如圖4.21：
x = seq(0,50, by = 0.1)
df<-data.table::data.table(x,
y = pexp(x, rate=1/15))
ggplot(df, aes(x=x, y=y)) +
geom_line(color="#4013ae")
Figure 4.21: 指數分配累積機率密度圖
- 圖4.22顯示x從0到30的累積機率：
x = seq(0,60, by = 0.1)
df<-data.table::data.table(x,
y = dexp(x, rate=1/15))
plot(df$y ~ df$x, type='l', xlab='x', ylab=expression(paste(f(x))))
cord.2x <- seq(0, 30, 0.01)
cord.2y <- dexp(cord.2x , 1/15)
polygon(c(0, cord.2x, 30), c(0, cord.2y, 0), border=NA,
col=col2rgb("yellow", 0.3))
Figure 4.22: 指數分配累積機率區域圖
- 我們可以用
R計算區域面積如下：
lambda=1/15
f <- function(x){
lambda*exp(-lambda*x)
}
integrate(f, 0, 30)0.8647 with absolute error < 9.6e-15
qexp 回傳機率密度對應的百分位，或者\(F^{-1}(X=x_{0})=\frac{-\rm{ln}(1-p)}{\lambda}\)由低而高如圖4.23：
p = seq(0, 1, by = 0.1)
df<-data.table::data.table(p,
y = qexp(p, rate=1/15))
ggplot(df, aes(x=p, y=y)) +
geom_line(color="#4013ae")
Figure 4.23: 指數分配百分位圖
- 最後，
rexp 回傳指數分配\(\lambda\)值的隨機亂數，可排列成為機率密度圖4.24：
set.seed(02139)
N=1000
df<-data.table::data.table(expo = rexp(N, rate=15))
ggplot(df, aes(x=expo)) +
geom_histogram(aes(y=..density..), color='black', fill="white",
binwidth = 0.005, stat='bin') +
geom_vline(aes(xintercept=mean(expo)),
color="red", linetype="dashed", size=1) +
geom_density(alpha=.2, fill="#FF6666")
Figure 4.24: 指數分配機率密度圖
- 如果計算\(\lambda\)的平均值，會發現大約是0.066，也就是\(\frac{1}{15}\approx 0.066\)。換句話說，平均每分鐘載到0.066個乘客，也就是平均每15分鐘載到一位乘客。
mean(df$expo)[1] 0.06533
☛假設某位外送員，平均每5分鐘接到訂單。請問在未來1小時內，接到至少3張訂單、至多8張訂單的機率有多高？平均5分鐘有一張訂單，也就是1分鐘有0.2張訂單，1小時則有12張訂單。
5 描述統計
以最有效率的方式描述人口或是特定群體的量化或類別變數的重要特徵
分析單位可以是個體或是群體。
群體的描述統計，例如：國家的人均所得、房價中位數、收入不均指數（吉尼係數）、生育率\(\ldots\)等等。
個體的描述統計，例如：個人的性別、教育程度、滿意度、收入、政治態度、血糖指數\(\ldots\)等等。
可分為中央趨勢以及離散程度兩個面向。
可以以文字或者是圖形表示。
例如
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
- 眾數
- 中位數
- 百分位數
- 平均數
5.1 眾數(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資料，找出
5.2 中位數(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 的關係。圖 5.1顯示，生師比越低、低社會地位人口比例越低，房價的中位數越高。
library(ggplot2)
ggplot(data=MASS::Boston, aes(y=medv, x=ptratio)) +
geom_point(aes(color=lstat))
Figure 5.1: 波士頓各區的房價中位數與生師比及低社會地位人口比例散佈圖
5.3 中位數計算方式
- 當數列的數目是奇數，中位數是第\(\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這筆資料中，學生的中位數成績是多少？
5.4 四分位數(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)
quantile(X, c(.25, .5, .75), type=1)25% 50% 75% 1 1001 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)
quantile(m, c(0.25, 0.5, 0.75), type = 1)25% 50% 75% 250 311 420
- 例如我們模擬兩筆資料，平均數相同，但是離散程度不同， 以同樣的方法計算25, 50, 75分位，會得到不同的25及75分位。這個例子顯示即使中位數相同，但是資料分佈不同，25及75分位就不同。
set.seed(1000)
tmp1<-rnorm(100, 40, 10)
quantile(as.integer(tmp1), c(0.25, 0.5, 0.75), type=1)25% 50% 75% 34 40 45
set.seed(1000)
tmp2<-rnorm(100, 40, 14)
quantile(as.integer(tmp2), c(0.25, 0.5, 0.75), type=1)25% 50% 75% 31 40 47
5.5 百分位數(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
quantile(dt$Score, c(0.25, .75, 0.9), type=1)25% 75% 90% 75 88 92
quantile(dt$Score, c(0.25, .75, 0.9), type=4)25% 75% 90% 72.5 88.0 91.4
5.6 比較SPSS與R的輸出
R的輸出跟SPSS類似，我們可以加以對照（圖5.2）。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)
quantile(scores, c(.25,.5,.75,.9), type=6)25% 50% 75% 90% 41.75 64.00 71.75 78.00
v1<-here::here('Fig','v1_quantile.png')
knitr::include_graphics(v1)
Figure 5.2: 四分位統計
5.7 平均數
- 平均數衡量資料的中心位置，可以想成是觀察值的平衡點：比平均值大的數的總和等於比平均值小的數的總和的絕對值。
- 平均數會受到極端值的影響，可以用
trim刪除若干百分比的數。 - 可分為算術平均數跟加權平均數。
5.7.1 算數平均數：
- 公式如下：\[\bar{y}=\frac{\sum y}{n}\]
- \[y={6, 7, 8, 8, 9, 10, 13, 15, 16, 45}\]
- \[\bar{y}=\frac{\sum (6+7+\cdots , +45)}{10}=13.7\]
5.8 加權平均數
- 在不知道個別觀察值，只知道分組的個案數跟平均數，我們可以假設觀察值分為\(k=1\cdots k\)個組，每一組有\(y_{1}\),\(y_{2}\),\(\ldots\) 人，每一組平均數為\[\bar{y_{1}}, \bar{y_{2}},\cdots\], 則全體的平均數為： \[\bar{y}=\frac{\sum y_{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))
A$station1 [1] 25 33 44
$station2 [1] 43 66 78 81
$station3 [1] 90 76 105 110 121
groupn=sapply(A, length); groupnstation1 station2 station3 3 4 5
submean=sapply(A, mean); cat("air pollution of each station:", submean,"\n")air pollution of each station: 34 67 100.4
cat("Sum of air pollution=", sum(groupn*submean),"\n")Sum of air pollution= 872
totaln=sum(sapply(A, length));
cat("Total number of stations=",totaln,"\n")Total number of stations= 12
totalair=sum(sapply(A, sum));
cat("Sum of air pollution=",totalair,"\n")Sum of air pollution= 872
cat("average air pollution=", totalair/totaln)average air pollution= 72.67
sapply()函數可套用函數在列表的每一個向量。 以分組平均數計算總平均數：
#計算三個組各組平均數
mn=sapply(A, mean); mnstation1 station2 station3 34.0 67.0 100.4
#平均數乘以每一組的個案數
mn*groupnstation1 station2 station3 102 268 502
#總和除以全部個案數
sum(mn*groupn)/totaln[1] 72.67
5.9 偏態(skewness)
- 正偏：右邊的尾巴較左邊長，眾數小於中位數而平均數大於中位數，偏態係數大於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 5.3: 偏態統計
5.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)
skewness(hsb2$write, type=1)[1] -0.4784
skewness(hsb2$write, type=2)[1] -0.482
- 首先是Stata偏態計算結果，如圖5.4：
stata<-here::here('Fig','write_stata.png')
knitr::include_graphics(stata)
Figure 5.4: Stata的偏態統計
- 再來是SPSS偏態計算結果，如圖5.5：
spss<-here::here('Fig','write_spss.png')
knitr::include_graphics(spss)
Figure 5.5: SPSS的偏態統計
5.10 峰度 (kurtosis)
- 測量資料的分佈是高聳或是平坦
- 越平坦則兩邊尾部越長，越高聳則是靠近平均值的部份越集中
- 計算峰度的公式： \[\frac{m^4}{m^2_{2}}-3\] \[s2=\sum (x_{i}-\bar{x})^2\] \[s4=\sum (x_{i}-\bar{x})^4\] \[m2=\frac{s2}{n}\] \[m4=\frac{s4}{n}\]
- 不同統計軟體計算峰度的公式略有不同，如果用
R，可以用e1071這個套件裡面的kurtosis，用type指令選擇2可得到跟SPSS一樣的答案。 - Stata的計算峰度公式為 \[\frac{m^4}{m^2_{2}}\]
- 樣本數目越大，理論上各種計算方式的結果越接近。
5.10.1 使用語法計算峰度
s4=sum((hsb2$write-mean(hsb2$write))^4)
s4[1] 3577773
m4=s4/200
s2=sum((hsb2$write-mean(hsb2$write))^2)
m2=s2/200
m4/m2^2[1] 2.239
5.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)[1] -0.7615
kurtosis(hsb2$write, type=2)[1] -0.7502
- 其次是Stata，如圖5.6：
stata<-here::here('Fig','write_stata.png')
knitr::include_graphics(stata)Figure 5.6: Stata的峰度統計
- 最後是SPSS，如圖5.7：
spss<-here::here('Fig','write_spss.png')
knitr::include_graphics(spss)Figure 5.7: SPSS的峰度統計
- 以上的結果顯示我們可以互相比較三種軟體計算的結果，
R的彈性比較大。
6 離散
6.1 範圍
- 範圍(range)：最大值及最小值的差距，也稱做全距。
- 若是常態分佈，範圍約等於六個標準差。
range(hsb2$write)[1] 31 67
range(hsb2$math)[1] 33 75
- 平均數相同，範圍可能不同；有的變數比較離散。
- 全距相同，但是離散程度可能不同，請見圖 6.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 6.1: 相同全距離散程度不同的機率密度圖
- 全距容易受到最大值與最小值的影響；最大值會拉大全距，離散程度變大。
6.2 四分位距
- 第一個四分位跟第三個四分位之間的差距。不受到極端值的影響。
- 可用在箱型圖或盒鬚圖 6.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 6.2: 學生成績盒鬚圖
- 可以看到，盒鬚圖的下緣約等於\(Q_{1}\)，上緣則是\(Q_{3}\)，中間的線則是中位數。盒鬚圖的中間部分是四分位距，也就是71.75-41.75=30。
☛請計算MASS::Animals這筆資料中的腦容量的四分位距。
6.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}\)
- 大於或等於0
- 如果樣本來自二元分佈，即0,1，則標準差為：
\[\sqrt{\frac{np(1-p)}{n-1}}\] - \(p\) 是事件發生的機率。
6.3.1 常態分佈的樣本標準差
如果隨機變數屬於常態分配，大部分的值應該聚集在平均值加減一個標準差的範圍內，因此，樣本標準差的大小特別重要。
當樣本來自於常態分配的母體，利用微積分可求出平均數的加減1個標準差包含約\(68\%\)的樣本。2個標準差包含約\(95\%\)的樣本。3個標準差包含約\(99\%\)的樣本。
先寫程式計算標準差，再用
R的sd() 函數驗證：
#hsb2 data
#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}}\)。
6.4 標準差的特性
- 改變樣本的單位，標準差也會改變
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})\]
7 標準常態分佈的Z值
\[f(Z)=\frac{1}{2\pi}e^{-\frac{Z^2}{2}}\quad Z=\frac{X-\mu}{\sigma}\]
- 標準常態分佈的性質為：
\[Z\sim N(0,1)\]
7.1 Z值或分數
- \(Z\)值可幫我們瞭解觀察值在資料中的相對位置。計算\(Z\)的公式為：
\[Z=\frac{x-\bar{x}}{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)-pnorm(-3,0,1)[1] 0.02365 - 可畫圖如圖7.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 7.1: 標準常態分佈曲線下的2.5%區域
★請問在
- 要求\(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
- 可畫圖如圖7.2：
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 7.2: 標準常態分佈曲線下的特定區域
★ 有一位員工的今年月薪為8.5萬，去年則為8萬。今年的全體薪水標準差為2.3萬，平均值為6.4萬，而去年的全體員工薪水標準差為2萬，平均值為6.2萬。請問該員工月薪相較於全體員工有增加嗎？
\[z_{1}=\frac{8.5-6.4}{2.3}=0.91\] \[z_{2}=\frac{8-6.2}{2}=0.75\] - 因為\(z_{1}\geq z_{2}\)，因此該員工月薪相較於全體員工有增加。
8 作業
8.1 隨機變數作業
2.1 \(f(x)=\frac{x^2-1}{6}\quad x=-2,1,2\) 2.2 \(f(x)=\frac{(x-1)^2}{15}\quad x=-2,-1, 0, 1,2\)
R的integrate函數）
8.2 描述統計作業
score中位數、90百分位數以及男性跟女性的平均數：
score平均數以及標準差：
UsingR套件中的faithful資料，請問要看噴泉最少要等幾分鐘？平均要等幾分鐘？最多跟最少等的時間差距多少分鐘？
airquality這筆資料的Wind這個變數，計算前後兩個資料點的差異，以分析兩天之間風速的差異。例如：
A <- c(35, 61, 69)
d.A <- c(26, 8)airquality這筆資料的Wind的偏態為何？峰度是多少？
councilor這筆資料中，請問平均工程預算是多少？樣本標準差多少？
2008Election這筆資料，請問馬英九的得票數的25分位數、中位數、75分位數分別是多少？請問25與75分位數之間差別多少？
9 更新內容日期
today <- Sys.Date()
today <- format(today, '%m/%d/%Y')
cat('最後更新日期', today )最後更新日期 04/15/2022