社會科學統計方法
隨機變數、機率原理
隨機變數、機率原理
蔡佳泓
2024
1 樣本空間
- \(S\)代表樣本空間(sample space)，表示所有可能結果(outcome)的集合。 例如，從沒有鬼牌的52張撲克牌中隨機抽出一張，一個可能的樣本空間是數字（A到K）（包括13個元素）。
- 有些實驗有兩個或多個可能的樣本空間。例如，從沒有鬼牌的52張撲克牌中隨機抽出一張，一個可能的樣本空間是數字（A到K），另外一個可能的樣本空間是花色（黑桃，紅桃，梅花，方塊）（包括4個元素）。
- 事件(event)代表單一、超過一個、或者完全沒有結果的集合。樣本空間的任何一個子集稱為事件。
- 例如屬於哪一個年級是一個事件，從1到5可能發生。
- 也可以是兩個複合的情形，例如又是1年級又是商學院的事件。
1.1 集合
集合類似事件，可以代表事件或者結果。例如6都的集合寫成:\(\{新北市、台北市,\ldots, 高雄市\}\)，太陽系行星的集合:\(\{\text{Mars},\text{Venus},\ldots,\text{Naptune}\}\)
集合論(set theory)有許多定義，例如A, B是兩個互斥(mutually exclusive)集合，代表這兩個集合之間沒有共通的事件。例如22個縣市跟太陽系行星沒有共通事件。如果C是非6都的縣市，那麼\(A\cup C\)等於22個縣市的樣本空間。但是A與C是互斥集合。
樣本空間可以是有限的集合，例如一個骰子所有出現的點數、一年之中所有的月份、寫成：
\[ S=\{\omega, 1, 2, 3, \ldots\} \]
- 樣本空間也可以是兩個樣本空間的組合，例如擲兩顆骰子得到的所有可能的點數，寫成：
\[ S=D\times D=\{(x,y)\hspace{0.5em}|\hspace{0.5em}x\in D, \hspace{0.5em}y\in D\} \]
其中\(D\)本身就是一個樣本空間。
樣本空間可以是無限樣本點的集合，例如打者的打擊率，寫成：
\[ S=\{x\in R\hspace{0.5em}|\hspace{0.5em}0\le x\le 1\} \]
- 常用Venn圖表示樣本空間與事件的關係，如圖1.1：
library(limma)
g <- cbind(
A = c(F, T, T, F, T, F, F, T),
B = c(F, F, F, T, T, F, F, T))
d <- vennCounts(g)
vennDiagram(d, circle.col = c('#0eccc1', '#aae010'))
Figure 1.1: Venn圖
在這個圖中，A是\(\{F, T, T, F, T, F, F, T\}\)，B是\(\{F, F, F, T, T, F, F, T\}\)，所以交集是2個\(\{T,T\}\)，至少有一個\(\{T\}\)是5個，都沒有是3個\(\{F,F\}\)。整個樣本空間有8個樣本。
可以進一步讀資料畫圖，也可以結合ggplot2以及vennCounts。
例如抽出一張撲克牌，有四種花色的樣本空間：\(S=\{\)紅心, 鑽石, 黑桃, 梅花\(\}\)。如果是不同花色與點數的樣本空間，總共有52個元素：\(S=\{\)紅心1, 紅心2,, 梅花K\(\}\)。n(S)=52。
自然萬物可以分為有機物與無機物，碳與氫組成的是有機物，例如動、植物是有機物，石頭、水是無機物，所以自然的樣本空間可以表示為：\(S=\{有機物、無機物\}\)。n(S)=2。
樣本空間內每一種可能的結果稱為元素或樣本點。
事件 (event) 為特定結果所形成的集合，事件為 \(S\) 的子集,以大寫字母 (\(A, B, C, \ldots\)) 表示，可分為以下兩種類型:
事件中只有包含一個元素稱作簡單事件 (simple event, 亦稱為樣本點) 例如投擲一顆骰子得點數 3。
事件中包含兩個以上的元素稱作複合事件 (compound event)，例如投擲兩顆骰子得點數和為 3 或 5 或 7。
若 A, B 兩事件滿足 \(A\cap B=\emptyset\)時, 則稱此A事件與B事件為互斥事件 (disjoint events)。A, B兩個事件不可能同時發生，這兩個事件中的元素不會互相重疊。
例如學生的年級成為一個樣本空間：\(S=\{\)大一，大二，大三，大四\(\}\)。一名學生不可能同時是大一以及大二的學生，所以屬於大一以及屬於大二是互斥事件。但是一個學生可以同時主修政治學以及會計學，所以主修可能不是互斥事件。
1.1.1 實例
- 擲一枚硬幣，樣本空間\(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\)
1.2 計算
- 有序的組合稱為permutation，例如從A到E抽出3個字排列，不同順序視為不同的結果，有\(P(5,3)=\frac{5!}{(5-3)!}=60\)種排列方式。得到A, B, C這個排列方式的機率是\(\frac{1}{60}\)。
- 無序的組合稱為combination，例如例如從A到D抽出3個字排列，不同順序視為相同，\(C(5,3)=\frac{5!}{3!(5-3)!}=10\)。得到A, B, C、C, B, A、A, C, B、B, A, C、C, A, B其中一個排列方式的機率是\(\frac{1}{10}\)。
1.2.1 範例
假設有10個家庭參加旅行團，而旅館的某一層有10間房間排成一排，請問其中兩個家庭剛好被安排住在一起的機率多高？
10個房間隨機安排給10個家庭，所以有10階乘的排列方式：\(10\times 9\times,...,1\)。
其中2個家庭想住在一起，10個房間代表有9階乘的排列方式：\(9\times 8\times,...,1\)，也就是有\(9!\)的安排方式。
所以兩個家庭住在一起的機率是\(\frac{9!}{10!}\)。
因為兩個家庭可以住一左一右，所以再乘以2。\(2\times \frac{9!}{10!}=\frac{2}{10}=\frac{1}{5}\)。
1.3 獨立事件 (Independent Events)
如果事件 \(A\) 發生時與事件 \(B\) 沒有關係，也就是 \(P(A,B)=P(A)\cdot P(B)\), \(\Pr(A|B)=\Pr(A)\)，這兩個事件稱為獨立事件。
例如：擲一顆骰子兩次，第一次得到6的機率與第二次得到6的機率相互獨立。兩次都得到6的機率是\(\frac{1}{6}\cdot\frac{1}{6}=\frac{1}{36}\)
抽一張撲克牌，放回去之後再抽一次，第一次抽到紅心的機率與第二次抽到紅心的機率相互獨立。但是如果抽出不放回，那麼這兩次抽到紅心的機率並不獨立，第二次抽到紅心的機率受到第一次抽到的結果影響。
問題：班上有20位同學，請問兩個同學同月同日生的機會多高？
因為「不同」生日是獨立事件，假設第一個同學生日是1月1日，第二位同學跟她不同生日，也許是1月2日，所以計算連續20位同學不同生日的機率為：
\[ P=\frac{365}{365}\times \frac{364}{365}\times \frac{363}{365}\times \ldots \times \frac{346}{365}=\frac{365\times 364\times \ldots \times 346}{365^{20}} \]
# Create a sequence of numbers from 365 to 346
numbers <- 365:346
# Compute the product of the numbers
result <- prod(numbers)
cat('P=', result/365^20)P= 0.5886
- 如果P代表沒有人同一天生日，那麼1-P等於有人同一天生日，所以：
\[ 1-P=0.414 \]
20人當中，出現同一天生日的機率高達4成。當然，人數越多，同一天生日的機率越高。所以遇到同一天生日的人其實不用太意外。詳細過程可參考這篇文章。
某個神射手的三分球命中率是0.56，請問他出手6次，前3次都沒投進，後3次連續投進的機率多高？
\[ (1-0.56)^3*0.56^3\approx 0.014 \]
- 兩個事件如果互相獨立，代表知道一個事件無法解釋另一個事件。例如台灣今天下雨跟明天英國下雨應該是互相獨立。高雄市路上的車輛跟台北市西餐廳的人數應該也是互相獨立。教育程度跟有沒有擁有汽車駕照可能是互相獨立，因為駕駛人分佈在各個社會階層。但是有沒有擁有護照可能不是互相獨立，因為教育程度越高，出國旅遊、工作的機會可能越多。
1.4 獨立事件實例
從52張撲克牌抽到紅色的牌與抽到有臉的牌(J, Q, K)是否為獨立事件?
如果是獨立事件，\(\Pr(A,B)=\Pr(A)*\Pr(B)\)
- 已知抽到紅色牌的機率\(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.5 互斥事件 (Mutually Exclusive/Disjoint Events)
- 兩個事件永遠不可能同時發生，稱為互斥事件。用Venn圖表示兩個事件沒有重疊，如圖1.2。
library(eulerr)
set.A <- c(LETTERS[1:6])
set.B <- c(LETTERS[7:12])
intersect(set.A, set.B)character(0)
plot(euler(list(A = set.A, B = set.B)),
fills = c("#aabb22", "#eeee00"),
labels = c("Set A", "Set B"))
Figure 1.2: 互斥事件
如果不是互斥事件，兩個事件的聯集是兩個事件的機率相加扣掉交集。
非互斥事件表示如圖1.3：
library(eulerr)
set.A <- c(LETTERS[1:6])
set.B <- c(LETTERS[1:3], LETTERS[7:9])
intersect(set.A, set.B)[1] “A” “B” “C”
plot(euler(list(A = set.A, B = set.B)),
fills = c("#11cccc", "#11ee00"),
labels = c("Set A", "Set B"))
Figure 1.3: 聯集事件
- 例如擲一個6面的骰子，A事件是出現偶數點數，\(A:\{2,4,6\}\)，B是至少4點，\(B:\{4,5,6\}\)，這兩個事件有交集，所以發生兩個事件的機率要扣掉交集：
\[ P(A\cup B)=P(A)+P(B)-P(A\cap B)=\frac{3}{6}+\frac{3}{6}-\frac{2}{6}=\frac{2}{3} \]
- 兩個事件永遠不會同時發生，聯合機率為0。兩個互斥事件的發生機率 \(P(A\cup B)=P(A)+P(B)\)。
- 例如，一個學生不可能同時為大一學生，也是大四學生，所以是互斥事件。但是一個學生可以同時為大一學生也是男生。
- 如果全部學生有1000人，可以分為大一、二、三、四，每個年級有同樣多的學生，隨機抽一個學生，該學生是大一生的機率是0.25，該學生是大一生或者大四生的機率是0.25+0.25=0.5。
- 例如，一個生物可能是動物，也可能是植物，但不可能同時是動物跟植物，例如珊瑚是動物，不是植物。
- 例如一個事件是奇數另一個事件是偶數: \(C=\{1,3,5\}, D=\{2,4,6\}\). \(C\) 與 \(D\) 事件為互斥。 \(C\cap B=\emptyset\). \(P(C,D)=0\)
- 互斥事件是永遠不會同時發生的事件，獨立事件則是不相關的事件，但是有可能同時發生。
- 例如，企鵝有很多品種，皇帝企鵝、巴布亞企鵝等等，每一個品種都有雄性、雌性，而且分佈平均，沒有那一種企鵝只有雄性，另一種企鵝都是雌性。因此，企鵝品種跟企鵝的性別是獨立事件。
- 兩個互斥事件如果發生其中一個，另一個就不發生，這兩個事件稱為互補事件(complement)。例如一個硬幣出現正面，就不會發生出現反面的事件。正面與反面是擲硬幣的樣本空間。
2 機率
給每一個事件一個數字\(P(A)\)，\(P:S\rightarrow R\)，其中，
- \(P(A)\geq 0\quad \forall A\subset S\)
- \(P(S)=1\)
- 對不交集的事件\(A_{1},A_{2},\ldots\)，\(P(\cup A_{i})=\sum P(A_{i})\)。例如國務院有東亞所、俄羅斯所、外交系、外交所等等9個系所與學程，如果運動會抽一個單位代表國務院，東亞所或者俄羅斯所被抽到的機會是\(\frac{1}{9}+\frac{1}{9}=\frac{2}{9}\)
以大樂透的中獎機率為例，我們可以先計算所有號碼的組合，每一組號碼的中獎機率為： \[ P(A)=\frac{n(A)}{n(S)}=\frac{n(A)}{n} \]
其中A代表頭獎的中獎號碼，n(A)代表A這個集合裡面的元素的數目，例如有三個人買了同樣的號碼，n(A)=3。n(S)代表所有號碼的組合數，例如大樂透是六個01到49的號碼，六個全對得到頭獎，所以一共有\(49\times 48\times\ldots\times 44=10068347510\)組號碼。
機率可以定義為無數次實驗之後，某事件A出現的相對次數，\(P(A)=\frac{n(A)}{n(S)}\)。隨機實驗是一個過程，在實驗前我們知道所有可能出現的結果，但是不確定是否會出現，也不知道會出現幾次。一個實驗或是隨機變數中所有可能出現的結果 (outcomes) 所形成的集合稱為樣本空間\(S\)。
2.1 機率公設
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)\)
☛請問這兩個函數是否符合機率公設？ \[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\)。
2.2 聯合機率 (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位黑人大法官，也就是\(\Pr(\text{Black, Female})=1/116\approx 0.8\%\)。
我們統計本國民眾各個年齡層、教育程度以及性別的結婚人數，為了節省篇幅，只看博士學位的男性民眾部分。一共有829位具有博士學位的男性民眾結婚，全部具有博士學位的結婚民眾則有1078位。25至29歲博士佔男性博士且結婚的民眾\(4.58\%\)，也就是P(25_29, 男, 博士)=\(4.58\%\)，但是佔全部具有博士學位的結婚人數的\(3.52\%\)，以此類推。
file <- here::here('data', 'opendata109m113.csv')
dt <- read.csv(file, sep=',', header = TRUE)
library(dplyr); library(janitor)
dt <- dt[-1,]
tab1 <- dt %>%filter(nation=='本國',sex == '男', edu == '博士畢業') %>%
summarise(total1=sum(as.numeric(marry_count)))
total1 <- tab1[1,1]
tab2 <- dt %>% filter(nation=='本國') %>%
group_by(edu, age) %>%
summarise(n = sum(as.numeric(marry_count))) %>%
mutate(total2=sum(n)) %>%
mutate (percentage2 = 100*n/total2) |>
mutate (percentage1 = 100*n/total1)
knitr::kable(tab2[1:12, ], format = 'simple',
booktabs=T,caption="109年本國民眾年齡層以及教育程度的結婚統計")| edu | age | n | total2 | percentage2 | percentage1 |
|---|---|---|---|---|---|
| 博士畢業 | 15～19歲 | 0 | 1078 | 0.000 | 0.000 |
| 博士畢業 | 20～24歲 | 0 | 1078 | 0.000 | 0.000 |
| 博士畢業 | 25～29歲 | 38 | 1078 | 3.525 | 4.584 |
| 博士畢業 | 30～34歲 | 348 | 1078 | 32.282 | 41.978 |
| 博士畢業 | 35～39歲 | 380 | 1078 | 35.251 | 45.838 |
| 博士畢業 | 40～44歲 | 156 | 1078 | 14.471 | 18.818 |
| 博士畢業 | 45～49歲 | 56 | 1078 | 5.195 | 6.755 |
| 博士畢業 | 50～54歲 | 38 | 1078 | 3.525 | 4.584 |
| 博士畢業 | 55～59歲 | 23 | 1078 | 2.134 | 2.774 |
| 博士畢業 | 60～64歲 | 16 | 1078 | 1.484 | 1.930 |
| 博士畢業 | 65歲以上 | 23 | 1078 | 2.134 | 2.774 |
| 博士畢業 | 未滿15歲 | 0 | 1078 | 0.000 | 0.000 |
2.3 邊際機率（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')
flextable::proc_freq(councilor, "Year", "unit")Year | unit | ||||
|---|---|---|---|---|---|
公園處 | 新建工程處 | 水利處 | Total | ||
2015 | Count | 1 (10.0%) | 1 (10.0%) | 2 (20.0%) | |
Mar. pct (1) | 100.0% ; 50.0% | 100.0% ; 50.0% | |||
2016 | Count | 8 (80.0%) | 8 (80.0%) | ||
Mar. pct | 100.0% ; 100.0% | ||||
Total | Count | 1 (10.0%) | 8 (80.0%) | 1 (10.0%) | 10 (100.0%) |
(1) Columns and rows percentages | |||||
#tu<-table(councilor$Year, councilor$unit)
#tu
#margin.table(tu,1)/sum(tu)- 計算得知，2015年與2016年的邊際機率分別為0.2及0.8。三個單位則分別是0.1,0.8,0.1。
2.4 條件機率（Conditional Probability）
- 隨機變數Y出現y而另一個變數X出現x的機率，或者當事件B發生時事件A也發生，表示為\(P(A|B)\). \(A\) 發生是以\(B\)發生為前提。
- 條件機率可以用聯合機率以及邊際機率表示：
\[P(A|B)=\frac{P(A\cap B)}{P(B)}\]
- 以表2.2為例，學生而沒有債務違約的條件機率是95.7%(\(2817/2944\))，非學生而有債務違約的條件機率是97.1%。
library(dplyr); library(ISLR)
ftb <- flextable::proc_freq(Default, "default", "student")
ftb <- flextable::set_caption(ftb, "學生身份與債務違約交叉表")
ftbdefault | student | |||
|---|---|---|---|---|
No | Yes | Total | ||
No | Count | 6,850 (68.5%) | 2,817 (28.2%) | 9,667 (96.7%) |
Mar. pct (1) | 97.1% ; 70.9% | 95.7% ; 29.1% | ||
Yes | Count | 206 (2.1%) | 127 (1.3%) | 333 (3.3%) |
Mar. pct | 2.9% ; 61.9% | 4.3% ; 38.1% | ||
Total | Count | 7,056 (70.6%) | 2,944 (29.4%) | 10,000 (100.0%) |
(1) Columns and rows percentages | ||||
- 從上一個式子可以推導出：
\[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}\]
如果我們知道當A發生時，B發生的機率，而且知道A、B各自發生的機率，那麼根據貝氏定理，我們可以求出當B發生時，A發生的機率。
以表2.2為例，如果我們要計算\(\Pr(\text{default}=Yes\mid\text{student}=No)\)可以這樣計算：
\[\begin{align} \frac{\Pr(\text{student}=no\mid\text{defalut}=Yes)\Pr(\text{default}=Yes)}{\Pr({\text{student}=no)}} & = \frac{0.619\times 0.033}{0.706}\\ & \approx 0.029=2.9\% \end{align}\]
2.5 相依事件 (Dependent Events)
已知當兩個事件是獨立事件時，則： \[ P(A\cap B)=P(A)\cdot P(B) \]
又因為
\[ \frac{P(A\cap B)}{P(A)}=P(B\mid A) \]
- 所以如果兩個事件獨立，則：
\[ P(A\cap B)=P(A)\cdot P(B)=P(B\mid A)\cdot P(A)\\ \Rightarrow P(B)=P(B\mid A) \]
B的邊際機率等於條件機率。
例如在表2.2，\(\Pr(\text{student}=Yes,\text{defalut}=Yes)\neq \Pr(\text{student}=yes)\times \Pr(\text{default}=Yes)\)，所以這兩個事件是互賴事件。
例如在表2.3，\(P(\text{food}=rice)=P(\text{nature}=insects\mid \text{food}=rice)\)，所以這兩個事件是獨立事件。
library(flextable)
me <- data.frame(nature=c(rep("insects", 100), rep('mammals', 100)),
food=rep(c(rep("rice", 50), rep("noodle", 50)),2))
ftab <- flextable::proc_freq(me, "nature", "food")
ftab <- flextable::set_caption(ftab, "觀察自然現象與喜歡吃的主食")
ftabnature | food | |||
|---|---|---|---|---|
noodle | rice | Total | ||
insects | Count | 50 (25.0%) | 50 (25.0%) | 100 (50.0%) |
Mar. pct (1) | 50.0% ; 50.0% | 50.0% ; 50.0% | ||
mammals | Count | 50 (25.0%) | 50 (25.0%) | 100 (50.0%) |
Mar. pct | 50.0% ; 50.0% | 50.0% ; 50.0% | ||
Total | Count | 100 (50.0%) | 100 (50.0%) | 200 (100.0%) |
(1) Columns and rows percentages | ||||
- 例題：有10位學生搭綠1，15位學生搭棕18，如果有6位非學生搭棕18，請問當有幾位非學生搭綠1時，學生身份與搭的公車是獨立事件?
3 全機率定理（Law of Total Probability）
- 假設我們有兩個事件A, B，其中B所在的樣本空間可以分為許多事件，\(P(B_{1})+P(B_{2})+,\ldots ,+P(B_{n})=S\)（樣本空間），那麼以下的關係成立：
\[ P(A|B)=P(A|B_{1})\cdot P(B_{1})+P(A|B_{2})\cdot P(B_{2})+,\ldots ,+P(A|B_{n})\cdot P(B_{n}) \]
- 全機率定理的基礎是樣本空間如果可以劃分為\(E_{1},\ldots, E_{n}\)，那麼另一個事件A跟這個樣本空間的交集寫成：
\[ P(A) = P(A\cap E_{1})\cup P(A\cap E_{2})\cup,\ldots,P(A\cap E_{n}) \]
- 例如A表示男性這一個事件，B與B’或者\(B^{c}\)表示成年與否這兩個來自年齡樣本空間的事件，\(P(A)\)可以寫成：
\[ P(A)=P(A\cap B)\cup P(A\cap B) \]
- 可以再切割年齡為成年、青少年、兒童等等階段，\(P(A)\)可以表示為：
\[ P(A)=P(A\cap B_{1})\cup P(A\cap B_{2})\cup,\ldots,P(A\cap B_{j})=\sum_{j=1}^nP(B_{j}\cap A) \]
要強調的是B的樣本空間應該可以切割成互斥且互補的事件，也就是窮盡所有可能出現的結果。
全機率定理是統計的重要原理，被廣泛應用在貝氏統計，以及機器學習等領域。
假設網路上只有三個冷凍水餃的賣家A、B、C，分別佔了40%、30%、30%的市場。物流公司從這三個賣家在每個月一日收到訂單的機率分別是5%、10%、15%。請問物流公司在每個月一日收到冷凍水餃訂單的機率用\(P(D)\)表示，如何計算\(P(D)\)？
因為訂單有可能來自A、B、C其中一家，但是這三家的市佔率以及接到訂單的機率都不同。以A家來說，在\(40\%\)市佔率的情況下，接到訂單機率是\(5\%\)，所以訂單來自A的機會是\(0.05\times 0.4\)。其他家也是如此。
所以物流公司收到訂單的機會應該是這三家公司的市佔率乘以接到訂單機率的總和，也就是：
\[\begin{align} P(D)=P(D|A)P(A)+P(D|B)P(B)+P(D|C)P(C)\\ & = 0.05\times 0.4+0.1\times 0.3+0.15\times 0.3 \\ & = 0.02+0.03+0.045 \\ & =0.095 \end{align}\]
- 另外一個例子跟醫學檢驗有關。假設有一台X光機，有90%的機率會正確驗出病人有結核病，有10%的機會沒有驗出來。如果病人並沒有結核病，有99%的機率正確地顯示病人沒有結核病，但是有1%的機會誤判有結核病。假設全台灣只有0.01%的人有結核病，隨機抽一位民眾去檢驗，請問機器顯示有結核病的機率多高？
4 隨機變數
4.1 定義
隨機變數的定義是一個轉換數字的函數，它的範圍是樣本空間。連續變數可以轉換任何數字，而不連續變數則轉換特定數字。
例如一個球員過去投進三分球的機率是45%，如果X代表他今天投n球所進的球，\(x=0\sim n\)，機率分別是多高？定義變數為\(X\sim P(x, n, p)\)，函數為：
\[ P(X=x) = \binom{n}{x}\times p^{x}\times (1-p)^{n-x} \]
- 當p=0.45, x=0, n=20，
\[ P(x=1)=\binom{10}{0}\times 0.45^0\times 0.55^{10}=0.002 \]
- 當p=0.45, x=1, n=20，
\[ P(x=1)=\binom{10}{1}\times 0.45^1\times 0.55^{10-1}=0.0207 \]
- 當p=0.45, x=2, n=20，
\[ P(x=2)=\binom{10}{2}\times 0.45^2\times 0.55^{10-2}=0.076 \]
- 依此類推，可以得到一個機率分佈：
# Define the parameters
n <- 10
p <- 0.45
# Generate the values of x
x_values <- 0:n
# Calculate the probabilities
probabilities <- dbinom(x_values, size = n, prob = p)
# Plot the distribution
plot(x_values, probabilities, type = "h", lwd = 2,
xlab = "Number of Successes (x)", ylab = "Probability",
main = "Binomial Distribution (n = 20, p = 0.45)")
當p=0.5，二元分佈的機率分佈將會是對稱。
更正式的寫法是對於任何實數x；
\[ f_{x}(X)=P(X=x)\hspace{1em}\forall\hspace{1em} x \in \mathbb{R} \]
- 與機率相同的隨機變數公設有：
- \(0\le f_{x}(X)\le 1\)
- \(\sum f_{x}(X_{i})=1\)
- \(P(X\subset A)=\sum_{A} f_{x}(X_{i})\)
5 間斷變數
5.1 間斷變數的期望值與變異數
5.1.1 期望值
以表5.1為例，計算丟擲兩枚硬幣出現正面的期望值。計算結果應為1。也就是說，我們預期出現正面的次數為1。
| 正面次數 | 相對次數或機率函數 | 相乘 |
|---|---|---|
| 0 | 0.25 | 0 |
| 1 | 0.5 | 0.5 |
| 2 | 0.25 | 0.5 |
| sum | 1 | 1 |
R模擬從兩個集合\(\{0,1\}\)各抽樣1000次，然後加總抽樣結果，並且畫成圖 5.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 5.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 5.2: 擲兩枚硬幣1000次的正面次數長條圖二
im1<-here::here('Fig','IMG_3167.JPG')
knitr::include_graphics(im1)
Figure 5.3: 打珠台
im2<-here::here('Fig','IMG_3168.JPG')
knitr::include_graphics(im2)
Figure 5.4: 打珠台事件與獎品
5.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})\)，就能計算變異數，以表5.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\)。我們可以求出二元分布的期望值\(E[X]=p\)，因為\(1\times p+0\times (1-p)=p\)。
- 變異數為\(p\times(1-p)\)。因為\(0\le p\le 1\)，所以變異數不會小於0。
★ 二項分佈(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- 可以畫圖表示如圖 5.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 5.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分佈是一段期間內發生事件的次數的機率分佈，例如從9點到10點進入郵局的人數，去年發生颱風的次數。
Poisson分佈的PMF寫成：
\[p(x=k;\lambda)=\frac{e^{-\lambda}\lambda^{k}}{k!}\quad \rm{for}\quad k=0,1,2,\ldots\]
- 其中Poisson分佈的平均值為\(\lambda\)，\(\lambda\)是在一定期間內的事件發生次數的平均數。圖5.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 5.6: 3種平均數的Poisson分佈
平均數\(\lambda\)越大，越接近常態分佈。
如果進行n次白努利試驗，結果會近似Poisson分佈。這是因為白努利分佈的期望值是\(n\times p\)，剛好是Poisson分佈的期望值\(\lambda\)。
set.seed(02138)
k = 1:20; p = 0.02
d = tibble(n = c(200, 400, 600, 1000, 1200, 1500))
random_binomial = function(k, n) {
rbinom(n, k, p)
}
d = d %>% mutate(randoms=purrr::map(n, random_binomial, k))
#d %>% purrr::unnest(randoms)
d %>% tidyr::unnest(randoms) %>%
ggplot(aes(x = randoms))+
geom_density(aes(fill = as.factor(n)))+
facet_wrap(~n, scales="free") +
theme(legend.position = 'none')
Figure 5.7: 同樣成功次數不同試驗次數的白努利分佈
- 當\(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.047165.1.3 Poisson分佈實例1
現代統計學 第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
- Poisson 分佈的累積機率密度函數如圖5.8：
# Define average rate (lambda)
lambda <- 10 # You can change this value
# Create vector of possible number of events (x) - adjust range as needed
x <- 0:30
# Calculate cumulative probabilities using ppois (Poisson CDF function)
y <- ppois(x, lambda = lambda)
#data frame
poissdf <- data.frame(x, y)
# Create the plot
p <- ggplot(aes(x = x, y = y), data = poissdf) +
geom_line(col = '#DD12BB') +
labs(x = "Number of Events",
y = "Cumulative Probability",
title = "Poisson Distribution CDF with lambda = 10") +
theme_bw()
p + ggplot2::geom_vline(xintercept = 10, linetype="dashed", color = "#DD1100")
Figure 5.8: Poisson 分佈的累積機率密度函數
- Poisson分佈的PDF是\(P(X=k, \lambda)\)，把k對應的機率累積，就可以得到圖 5.8 顯示的Poisson分佈的累積機率：\(f(X)=P(X>k|\lambda)\)。當\(k\)越大，\(P(X>k|\lambda)\)也越大，然後就接近1。曲線上的每一個點代表當X=k時，函數所累積的機率。
5.1.4 Poisson分佈實例2
現代統計學 第202頁，已知咖啡店平均每小時有50位顧客，請問每小時需要幾位服務人員，才不會讓每位服務人員需要服務超過4位客人的機率高於0.5？換句話說，在\(k>4\)的條件下，\(\lambda\)應該等於多少，才不會出現一個服務生必須服務超過4個客人？\(\lambda=50/b\)，\(b\)是服務生人數。圖 5.9顯示機率密度函數， 當k從0增加到20時，\(P(X=k\mid \lambda)\)在\(\lambda = 5\)時達到最高點，約0.175，而\(P(X=10\mid \lambda=5)\approx 0.018\)。也就是說，當有10位服務人員時，平均一位服務4人的機率為0.175。但是我們有興趣的是必須服務超過4人的機率有多高？
k<-seq(0, 20, by=1)
y <- dpois(k, lambda= 5)
plot(k, y, type='l', lwd=1.5, col='#1a00ab',
ylab=expression(paste('P(X=k|', lambda,')')))
Figure 5.9: Poisson機率分佈之機率密度圖
- 圖 5.10則顯示Poisson的累積機率分佈。
- 當\(P(X=0,\lambda=5)+P(X=1,\lambda=5)+\ldots+P(X=4,\lambda=5)=0.440\)，也就是1 - \(P(X>4|\lambda = 5) \approx 0.559\)。簡單來說，當\(\lambda=5\)時，平均服務人數超過4人的機率有55.9%。
showtext::showtext_auto()
#k
k<-seq(0, 20, by=1)
#cumulative density
y <- ppois(k, lambda = 5)
#Plot
plot(k, y, type='l', lwd=1.5, col='#e22e11', yaxt="n",
xaxt="n",
ylab=expression(paste('P(X>k|', lambda,')')), xlab = '平均服務的客人數')
axis(side = 1, labels = c('1','3','5','7', '9','11', '13',
'15', '17', '19'), at = seq(1,19, by=2))
axis(side = 2, labels = c('0.1','0.3','0.5','0.7', '0.9'), at = seq(0.1, 0.9, by = 0.2))
abline(v = 5, lty = 2)
abline(h= 1-ppois(4, 5), lty = 2)
Figure 5.10: Poisson機率分佈之累積機率密度圖
- 我們可以應用
ppois(k=4, \(\lambda\)) ，求出當平均值為\(50/b\)、k>4的時候，b應該等於多少，才會讓顧客被服務的平均人數至少為4的機率不到0.5。這裡要注意的是，ppois(k=4, \(\lambda\)) 是從0累加到4，也就是顧客被服務的平均人數0一直到4，平均人數越多應該需要越多服務人員。但其實我們關心的是超過4之後的機率是否小於0.5，因此要用1來減去累積機率，或者直接要R計算從曲線另一邊累積的機率。圖 5.11 顯示當b = 10，\(P(X\le 4\mid 5)\)的累積機率小於0.5。當然，咖啡廳可以雇用更多服務人員，\(P(X>4\mid \lambda)\)的累積機率也越小，也就是不太可能出現一位服務人員需要服務超過4位客人的情況，但是可能會超出雇主的負擔。
showtext::showtext_auto()
b <- seq (4, 12, by = 1)
y<-c()
for (i in 1: 9){
y[i] <- 1-ppois(4, lambda=c(50/b[i]))
}
plot(b, y, type='l', lwd=1.5, col='#BB2200',
ylab=expression(paste('1-P(X>4|', lambda,')')), xlab = '服務生人數')
abline(h = 0.5, lty=2, col='blue')
abline(v = 10, lty=2, col='blue')
Figure 5.11: 超出特定平均值之累積機率密度圖
6 連續變數
身高、體重、體溫、收入、雨量、時間、距離等等是常見的連續變數，因為這些變數的個數是無限而且不可數。
連續變數定義為：如果有一個函數，對任何一個區間，也就是\(A\subset R\)，有如下的關係就是連續變數：
\[ P(X\subset A)=\int _{A}f_{x}(x)dx \]
\(f_{x}(x)\)是PDF。
連續變數的公設有：
- \(0\le f_{x}(X)\)
- \(\int f_{x}(X_{i})=1\)
- \(P(A)=P(a\le X\le b)=\int_{A} f_{x}(x)dx\)
連續變數常用的機率分配有：常態分佈、標準常態分佈、均等分佈、指數分佈等等。
定義：連續隨機變數\(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個值，然後畫長條圖如圖 6.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 6.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\]
- 例如常態分佈的部分面積如圖 6.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 6.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\]
- 所以連續變數的CDF有如下的性質：
- \(0\le F(x)\le 1\)
- \(\lim_{x\rightarrow\infty}F(x)=1\)
- \(\lim_{x\rightarrow -\infty}F(x)=0\)
- 往右連續（如果是不連續變數則有跳躍）
- CDF與PDF的關係為：
\[ F_{x}(x)=P(X\le x)=\int_{-\infty}^{\infty}f_{x}(x)dx\\ FX(x)=\frac{dF(x)}{dx}=f_{x}(x)^{-\infty} \]
- 期望值：
\[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}\]
6.1 連續變數的期望值與變異數
6.1.1 期望值
\[\boxed{E(X)=\int_{a}^{b}xf(x)dx=\mu} \]
- 在間斷變數的例子，變數有個別的變量，例如變數\(X\)有四個可能的值1, 2, 3, 4，發生的機率如表 6.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 6.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*}\]
- 圖 6.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 6.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 6.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計算。
6.1.2 變異數
- 在了解連續變數的平均數之後，我們想要知道變數的離散程度。連續變數的變異數的計算方式為：
\[\begin{align} V(X)=\int_{a}^{b}(x-\mu)^2\cdot f(x)dx \tag{6.1} \end{align}\]
- 方程式 (6.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{6.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)\)的期望值以及變異數？
6.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 6.6: 常態機率分佈
- 如果資料近似鐘形的常態分佈，可以在長條圖上面加上常態分佈曲線，如圖 6.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 6.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 6.8: 累積機率密度圖
R提供的函數畫累積密度曲線，如圖 6.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 6.9: R函數之累積密度曲線圖
- 圖 6.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 6.10: 機率分佈曲線下的一半區域
6.2.1 常態分配相關函式
R的常態分佈函式有dnorm 產生機率密度函數，也就是代入X之後的值。我們以IQ這個連續變數為例，最小為70、最高為150，平均值110、標準差則是任意的15。我們可以產生對應的機率密度，並且畫成圖6.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 6.11: IQ機率密度圖
pnorm 產生累積機率密度函數，也就是代入X之後產生從最小值一直累積到X的機率。畫成圖 6.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 6.12: IQ累積機率密度圖
pnorm 產生累積機率密度函數，qnorm 則是pnorm 的相反，也就是產生從最小值一直累積到X的機率所對應的值。畫成圖6.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 6.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 產生來自常態分佈的隨機亂數，我們產生三個不同樣本規模但是相同的平均數以及變異數的樣本，然後畫相對次數的長條圖 6.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 6.14: 三種樣本規模的長條圖
6.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{6.3} \end{equation*}\]
- 在方程式(6.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)，\)圖 6.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 6.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
- 圖 6.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 6.16: 一致機率分佈曲線下的區域
6.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}(b-a)^3 \tag{6.4} \end{align*}\]
\[\begin{equation*} f(x) = \begin{cases} \frac{1}{X} & 0\le x\le 1\\ 0 & \text{otherwise} \end{cases}\tag{6.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\]
- 或者根據方程式(6.4)：
\[\begin{align*} E(Z^2)-E(z)^2 & =\frac{1}{1-0}\frac{1}{3}(1-0)^3-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
6.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 6.17: 均勻分佈的累積機率圖
qunif 回傳均勻分佈的百分位如圖6.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 6.18: 均勻分佈的百分位圖
runif 產生均等分配的隨機變數，例如抽出位於(0,10)的1000個數字，繪成長條圖6.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 6.19: 均等分配隨機變數的長條圖
punif(1, 0, 10)[1] 0.1
6.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}}\]
- 累積機率密度函數則是：
\[F(X;\lambda)=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\]
- 在Poisson分佈中，第k次隨機事件與第k+1次隨機事件出現的時間間隔服從指數分布，長度為t的時間段內有一次隨機事件出現的機率等於
\[ \frac{e^{-\lambda t}(\lambda t)^1}{1!} = e^{-\lambda t}(\lambda t) \]
- 第k次隨機事件之後長度為t的時間段內，第\(k+n\)次 (\(n=1, 2, 3,\ldots\))隨機事件出現的機率等於\(1-e^{-\lambda t}\)。
6.4.1 指數分配期望值與變異數
\[E(X)=\frac{1}{\lambda}\]
\[V(X)=\frac{1}{\lambda^2}\]
6.4.2 指數分配的相關指令
dexp 回傳機率密度，以上一題為例，可假設\(0\leq X\leq 50\)，得到對應的機率密度由高而低如圖6.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 6.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\)，得到對應的機率密度由高而低如圖6.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 6.21: 指數分配累積機率密度圖
- 圖6.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 6.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}\)由低而高如圖6.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 6.23: 指數分配百分位圖
- 最後，
rexp 回傳指數分配\(\lambda\)值的隨機亂數，可排列成為機率密度圖6.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 6.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張訂單。
7 機率作業
X <- c(0, 1, 2, 3)
f <- c('1/2', '1/4', '1/8', '1/8')
df <- data.frame(X, f)
colnames(df) <- c('X', 'f(x)')
kableExtra::kable(df, format='simple')| X | f(x) |
|---|---|
| 0 | 1/2 |
| 1 | 1/4 |
| 2 | 1/8 |
| 3 | 1/8 |
\[f(x)=\frac{x^2-1}{6}\quad x=-2,1,2 \]
\[f(x)=\frac{(x-1)^2}{15}\quad x=-2,-1, 0, 1,2 \]
- 1. 擲3次硬幣，每次出現朝正面的事件。
- 2. 每次搭飛機遇到亂流的事件。
- 3. 3個人同時抽籤決定其中一人去倒垃圾時，每一個人抽到籤的事件。
- 4. 擲六面的骰子，第一次出現偶數的事件為A，第二次出現5或6的事件為B，這兩個事件是否為獨立事件？8 隨機變數作業
library(ggplot2)
uniform_Plot <- function(a, b){
xvals <- data.frame(x = c(a, b)) #Range for x-values
ggplot(data.frame(x = xvals), aes(x = x)) + xlim(c(a, b)) + ylim(0, 1/(b - a)) +
stat_function(fun = dunif, args = list(min = a, max = b), geom = "area",
fill = "green", alpha = 0.35) +
stat_function(fun = dunif, args = list(min = a, max = b)) +
labs(x = "\n u",
# title = paste0("Uniform Distribution \n With Min = ", a, " & Max = ", #b, " \n"),
y = "f(u) \n") +
theme(plot.title = element_text(hjust = 0.5),
axis.title.x = element_text(face="bold", colour="blue", size = 12),
axis.title.y = element_text(face="bold", colour="blue", size = 12)) +
geom_vline(xintercept = a, linetype = "dashed", colour = "red") +
geom_vline(xintercept = b, linetype = "dashed", colour = "red")
}
uniform_Plot(0, 10)
#Ans: f(x) = 1/10 or 0- 請列出此樣本空間。
- 如果隨機變數X表示加起來的點數，請列出每一個結果所對應的值。</h7>X <- c(20, 30, 40, 50, 60, 70, 80, 90)
f <- c(0.05,0.1,0.1,0.1,0.25,0.25,0.1,0.05)
df <- data.frame(X, f)
colnames(df) <- c('人數X', 'f(x)')
kableExtra::kable(df, format='simple')| 人數X | f(x) |
|---|---|
| 20 | 0.05 |
| 30 | 0.10 |
| 40 | 0.10 |
| 50 | 0.10 |
| 60 | 0.25 |
| 70 | 0.25 |
| 80 | 0.10 |
| 90 | 0.05 |
- 請問每一趟載客人數的期望值以及變異數為何？
- 如果每位乘客投幣15元，請問公車收到車資的期望值以及變異數為何？
R的integrate函數）
\[\begin{equation} f(x)= \begin{cases} x+0.5\quad -1<x<1\\ 0\quad \text{otherwise} \end{cases} \end{equation}\]
\[\begin{equation} f(x)= \begin{cases} 2x\quad 0<x<1\\ 0\quad \text{otherwise} \end{cases} \end{equation}\]
9 更新講義時間
library(lubridate)
# Get the current date and time
datetime <- Sys.time()
# Parse the current date and time using lubridate
current_datetime <- strftime(datetime)
# Print the parsed date and time
cat("最後更新時間:", current_datetime)最後更新時間: 2024-04-26 12:34:26