社會科學統計方法
隨機變數、機率原理
隨機變數、機率原理
蔡佳泓
4/18/2025
1 機率、事件、隨機變數
\[\begin{equation*} Y_{i}=x_{i}\beta+\epsilon_{i}\\ \epsilon_{i}\sim f_{n}(e_{i}|0,\sigma^2) \tag{1.1} \end{equation*}\]
\[\begin{equation*} Y_{i}\sim f_{n}(y_{i}|\mu,\sigma^2)\\ \mu_{i}=x_{i}\beta \tag{1.2} \end{equation*}\]
- 方程式(1.2)意思是\(Y_{i}\)的平均數\(\mu\)是\(x_{i}\)的函數，也就是
\[E(Y)=\mu=X\beta\]
- 假設\(Y\)從\(\theta\)與\(\alpha\)構成的分布抽樣得到，寫成方程式(1.3) :
\[\begin{equation*} Y\sim f(y|\theta,\alpha) \\ \theta=g(X,\beta) \tag{1.3} \end{equation*}\]
- 方程式(1.3)可以解讀成\(Y\)會隨著\(X\)特徵而變化，稱之為系統成分(systematic component)。而觀察到的\(y_{i}\)則是隨機成分(stochastis component)。在線性迴歸模型，\(\theta\)代表\(E(Y)\equiv \mu\)，\(\alpha\)是\(Y\)的重要資訊，但是我們更關心\(\theta\)。
\[\begin{equation*} \text{Pr}(y|M) \tag{1.4} \end{equation*}\]
- 其中\(y\)是資料，\(M\)是模型，包含\(\theta\), \(\alpha\)等參數，以及隨機變數\(Y\)等等。\(M\)代表不確定的程度。
- \({Pr}(y|M)\)大於0小於1。沒有單位。除了0跟1，其他的機率都有不確定性。
2 樣本空間
\(S\)代表樣本空間(sample space)，表示所有可能結果(outcome)的集合。樣本空間內每一種可能的結果稱為元素或樣本點。
例如，從沒有鬼牌的52張撲克牌形成一個的樣本空間是52張數字加上花色的組合。
有些實驗有兩個或多個可能的樣本空間。例如，從沒有鬼牌的52張撲克牌中隨機抽出一張，一個可能的樣本空間是數字（A到K），另外一個可能的樣本空間是花色（黑桃，紅桃，梅花，方塊）（包括4個元素）。
六都的樣本空間寫成:\(\{新北市、台北市,\ldots, 高雄市\}\)
全球前10大經濟體的樣本空間則是:\(\{美國、中國、歐盟、日本、德國、印度、英國、法國、加拿大、俄羅斯、義大利\}\)
事件 (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=\{\)大一，大二，大三，大四\(\}\)。一名學生不可能同時是大一以及大二的學生，所以屬於大一以及屬於大二是互斥事件。但是一個學生可以同時主修政治學以及會計學，所以主修不是互斥事件。
- 例如小學有六個年級，一個年級是一個事件，而一、二年級稱為低年級，也是一個事件。
- 也可以是兩個複合的情形，例如又是五年級、導師又是男性的事件。
3 集合
3.1 定義
\[V=\{x|x\hspace{.2em}\text{是義大利文的陽性名詞}\}\]
- 或者
\[V=\{libro, treno, ragazzo,\ldots\}\]
- 第二種方式是先定義\(A\)是所有義大利文的名詞，然後表示陽性名詞的集合：
\[V=\{x\in A|x\hspace{.2em}\text{是義大利文的陽性名詞}\}\]
3.1.1 例題
- 請問\(B=\{x\in R|x^2=4\}\)，\(B\)可以寫成？
3.2 集合論的概念
空集合：沒有任何元素的集合，寫成\(\emptyset\)。
相等集合：兩個集合如果有一樣的元素，稱為相等集合。
子集合：\(A\subseteq B\)表示\(A\)是\(B\)的子集合，也就是\(A\)的所有元素同時也是\(B\)的元素。理論上空集合是所有集合的子集合。此外，如果一個集合有\(n\)個元素，子集合會有\(2^{n}\)個。
集合之間有交集或者聯集。交集： \[A\cap B=\{x|x\in A \hspace{.6em}\text{and}\hspace{.6em}x\in B\}\]
聯集： \[A\cup B=\{x|x\in A \hspace{.6em}\text{or}\hspace{.6em}x\in B\}\]
- 例如：\(\mathcal{U}\)是113位立委的集合，\(A\)代表區域立委的集合，\(B\)是不分區立委的集合。\(A\cap B=\emptyset\)。而\(A\cup B=\mathcal{U}\)。
- 例如：\(A\)表示有手機門號的民眾，\(B\)表示有室內電話門號的民眾， \(C=\{x\in S|x\in A\hspace{.6em}\text{and}\hspace{.6em}x\in B\}\)表示同時有手機以及室內電話門號的民眾。\(C'\)表示兩種電話門號都沒有的民眾。
互斥：A, B是兩個互斥(mutually exclusive)集合，代表這兩個集合之間沒有共通的事件。例如22個縣市跟太陽系行星沒有共通事件。如果A代表六都的集合，B代表非六都的縣市的集合，那麼\(A\cup B\)等於22個縣市的樣本空間。但是A與B是互斥集合。
餘集合：A的餘集合代表除了A集合之外所有的事件。寫成：\(A'=\{x\in \mathscr{U}|x\notin A \}\)。
3.3 樣本空間與集合
- 樣本空間可以是有限的集合，例如一個骰子所有出現的點數，寫成：
\[ S=\{1, 2, 3, 4, 5, 6\} \]
- 樣本空間也可以是兩個樣本空間的組合，例如擲兩顆骰子得到的所有可能的點數，寫成：
\[ S=D\times D=\{(x,y)\hspace{0.5em}|\hspace{0.5em}x\in D, \hspace{0.5em}y\in D\} \]
其中\(D\)是一個樣本空間。
樣本空間可以是無限樣本點或者連續變數的集合，例如棒球打者的打擊率最小是0最大是1，寫成：
\[ S=\{x\in R\hspace{0.5em}|\hspace{0.5em}0\le x\le 1\} \]
例如：一個手機號碼在一個月內打了幾通是不連續變數的集合，而在一個月通話時間是連續變數的集合，因為通話時間從0到特定時間之間有無限的值。
用圖3.1表示連續變數的樣本空間及隨機變數與機率：
Figure 3.1: 連續變數的樣本空間
- 常用Venn圖表示樣本空間與事件的關係，如圖3.2：
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 3.2: 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。
3.3.1 樣本空間與集合的實例
抽出一張撲克牌，有四種花色的樣本空間：\(S=\{\)紅心, 鑽石, 黑桃, 梅花\(\}\)。如果是不同花色與點數的樣本空間，總共有52個元素：\(S=\{\)紅心1, 紅心2,, 梅花K\(\}\)。n(S)=52。
自然萬物可以分為有機物與無機物，碳與氫組成的是有機物，例如動、植物是有機物，石頭、水是無機物，所以自然的樣本空間可以表示為：\(S=\{有機物、無機物\}\)。n(S)=2。
擲一枚硬幣，樣本空間\(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\)
3.4 計算
- 有序的組合稱為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}\)。
3.4.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}\)。
4 事件
4.1 獨立事件 (Independent Events)
用機率的方式表示兩個事件為獨立事件： \(P(A,B)=P(A)\cdot P(B)\)，或者\(P(A\cup B)=P(A)\cdot P(B)\)， \(\Pr(A|B)=\Pr(A)\)，。
例如：擲一顆骰子兩次，第一次得到6的機率與第二次得到6的機率相互獨立。兩次都得到6的機率是\(\frac{1}{6}\cdot\frac{1}{6}=\frac{1}{36}\)。
例如：全公司25人抽獎，抽到頭獎的機率是\(\frac{1}{25}\)，連續兩年都抽到頭獎的機率是\(\frac{1}{25}\times \frac{1}{25}\)。
抽一張撲克牌，放回去之後再抽一次，第一次抽到紅心的機率與第二次抽到紅心的機率相互獨立。但是如果抽出不放回，那麼這兩次抽到紅心的機率並不獨立，因為第二次抽到紅心的機率受到第一次抽到的結果影響。
問題：班上有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 \]
- 兩個事件如果互相獨立，代表知道一個事件無法解釋另一個事件。例如台灣今天下雨跟明天英國下雨應該是互相獨立。高雄市路上的車輛跟台北市西餐廳的人數應該也是互相獨立。教育程度跟有沒有擁有汽車駕照可能是互相獨立，因為駕駛人分佈在各個社會階層。但是有沒有擁有護照可能不是互相獨立，因為教育程度越高，出國旅遊、工作的機會可能越多。
4.2 獨立事件實例
從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}\)，所以抽到紅色的牌與抽到有臉的牌是獨立事件。
4.3 互斥事件 (Mutually Exclusive/Disjoint Events)
- 兩個事件永遠不可能同時發生，稱為互斥事件。\(P(A\cup B)=0\)。用Venn圖表示兩個事件沒有重疊，如圖4.1。
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 4.1: 互斥事件
例如，一個學生不可能同時為大一學生，也是大四學生，所以是互斥事件。但是一個學生可以同時為大一學生也是男生。
如果全部學生有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\)。
如果不是互斥事件，兩個事件的聯集是兩個事件的機率相加扣掉交集。
非互斥事件表示如圖4.2：
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 4.2: 聯集事件
在這23個國家之中，屬於CPTPP的會員國的機率為\(\frac{12}{23}\)，RCEP則是\(\frac{15}{23}\)，交集則是\(\frac{4}{23}\)。
例如擲一個6面的骰子，A事件是出現偶數點數，\(A:\{2,4,6\}\)，B是至少4點，\(B:\{4,5,6\}\)，這兩個事件有交集，發生兩個事件的聯集是\(\{2,4,5,6\}\)，聯集的機率要扣掉交集部分\(\{4,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} \]
換句話說，擲一個6面的骰子擲出2,4,5,6的機率是\(\frac{2}{3}\)。注意交集的部分是\(\frac{2}{6}\)，因為在樣本空間1到6點之中，只有4,6兩個事件屬於交集的部分。
互斥事件是永遠不會同時發生的事件，獨立事件則是不相關的事件，但是有可能同時發生。
例如，企鵝有很多品種，皇帝企鵝、巴布亞企鵝等等，每一個品種都有雄性、雌性，而且分佈平均，沒有那一種企鵝只有雄性，另一種企鵝都是雌性。因此，企鵝品種跟企鵝的性別是獨立事件。
抽一張撲克牌是紅色字樣的機率是\(\frac{1}{2}\)，抽出6點的機率是\(\frac{4}{52}=\frac{1}{13}\)。因為這兩個事件是獨立事件，所以又是紅色又是6點的交集部分是\(\frac{1}{2}\times \frac{1}{13}=\frac{1}{26}\)。換句話說，52張撲克牌，抽到紅色6點只有2個事件。
兩個互斥事件如果發生其中一個，另一個就不發生，這兩個事件稱為互補事件(complement)。例如一個硬幣出現正面，就不會發生出現反面的事件。正面與反面是擲硬幣的樣本空間。
4.4 互賴事件 (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的邊際機率等於條件機率。
例如在表5.2，\(\Pr(\text{student}=\text{Yes},\text{defalut}=\text{Yes})\neq \Pr(\text{student}=\Yes})\times \Pr(\text{default}= \text{Yes})\)，所以這兩個事件是互賴事件。
例如在表4.1，\(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時，學生身份與搭的公車是獨立事件?
| 公車 | 綠1 | 棕18 | 和 |
|---|---|---|---|
| 學生 | 10 | 15 | 25 |
| 非學生 | x | 6 | |
| 和 | 21 |
- 假設全部人數\(N\)，由以上題目可知\(N=x+31\)。
- 獨立事件的條件是\(P(學生\cap 綠1)=P(學生)\times P(綠1)\)。
- \(P(學生)=\frac{25}{x+31}\)，\(P(綠1)=\frac{10+x}{x+31}\)，\(P(學生\cap 綠1)=\frac{10}{x+31}\)。x=4。
5 機率
隨機變數就是透過特定方式，給每一個事件一個數字\(X(A)\)，再轉換為機率\(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\)。
5.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)\)
☛請問這兩個函數是否符合機率公設？
\[\begin{equation} f(x)=\frac{x-1}{6}\quad x=1,2,3,4 \end{equation}\]
\[\begin{equation} f(x)=\frac{x^2}{2}\quad x=-1,0,1 \end{equation}\]
☛提示：(1)\(f(x)\)是否符合\(0\leq f(x)\leq 1\)；(2)代入x到\(f(x)\)，求得\(\sum f(x)=1\)。
5.2 聯合機率 (Joint probability)
\[f(x,y)=P(X=x,Y=y)\]
例如: 一張撲克牌，同時為紅色且為4點的機率為\(\Pr(4, \text{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\%\)，以此類推。
| 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 |
5.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。
5.4 條件機率（Conditional Probability）
- 隨機變數Y出現y而另一個變數X出現x的機率，或者當事件B發生時事件A也發生，表示為\(P(A|B)\). \(A\) 發生是以\(B\)發生為前提。
- 條件機率可以用聯合機率以及邊際機率表示：
\[P(A|B)=\frac{P(A\cap B)}{P(B)}\]
- 以表5.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發生的機率。
以表5.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}\]
6 全機率定理（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%的人有結核病，隨機抽一位民眾去檢驗，請問機器顯示有結核病的機率多高？
7 隨機變數
7.1 定義
隨機變數的定義是一個轉換數字的函數\(X\)，它的範圍是樣本空間。連續變數可以轉換任何數字，而不連續變數則轉換特定數字。
例如一個球員過去投進三分球的機率是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 \]
- 依此類推，可以得到一個機率分佈如圖7.1：
Figure 7.1: 進球分佈圖
當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})\)
8 不連續變數
8.1 以PMF與PDF表示事件發生機率
\[ 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} \]
- 或者寫成：
\[P(X=x)=\frac{1}{6},\quad x=1,2,3,4,5,6\]
- 而常態分佈的函數可寫成如下：
\[ 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的平方、次方、立方…連起來，用圖8.1表示：
#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 8.1: 自然對數的底
可以看出\(e\)有緩慢上升的趨勢，與常態分佈有關。
假設成人的身高成標準常態分佈，也就是標準差等於1而且平均數等於0，我們可以計算平均值左右一個Z值之間的身高的發生機率如圖 8.2。Z值計算方式：
\[Z=\frac{X-\mu}{\sigma}\]
# 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 = "", x = "Z value", y = "Density") +
theme_minimal()
# Display the plot
print(p)
Figure 8.2: 標準常態分佈與左右各一個Z值
8.1.1 範例
假設全體160,00名政大學生身高成常態分佈，平均值是167，標準差是5公分。今天抽樣統計，得知樣本當中身高最矮是160，最高175，請問這套樣本的身高在全體政大學生的當中，出現機率有多高？
用圖8.3表示政大學生身高分佈，以及樣本所佔的機率。
# Mean and standard deviation
mean_height <- 167
std_dev <- 5
lbound <- 160
ubound <- 175
# Generate data
height <- rnorm(1000, 167, 5)
# 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 >= lbound & x <= ubound),
aes(ymax = y, ymin = 0), fill = "#C11100", alpha = 0.3) +
labs(title = "",
x = "身高", y = "機率密度")
hp
Figure 8.3: 政大學生身高分佈
- 要得到落在中間的機率，用\(\textbf{pnorm}\)分別計算最小值到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 <- 167
std_dev <- 5
# 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.8644- 出現機率是0.86，也就是說這套樣本約代表政大86%的學生。圖8.3顯示樣本落在母體的中間偏左。
- 如果換成政大學生智商的分佈，我們可以假設成常態分佈嗎？
8.2 間斷變數的期望值與變異數
8.2.1 期望值
\[\begin{equation*} \boxed{E(X)=\sum_{i=1}^n x_{i}f(x_{i})=\mu} \tag{8.1} \end{equation*}\]
以表8.1為例，計算丟擲兩枚硬幣出現正面的期望值。計算結果應為1。也就是說，我們預期出現正面的次數為1。
| 正面次數 | 相對次數或機率函數 | 相乘 |
|---|---|---|
| 0 | 0.25 | 0 |
| 1 | 0.5 | 0.5 |
| 2 | 0.25 | 0.5 |
| sum | 1 | 1 |
R模擬從兩個集合\(\{0,1\}\)各抽樣1000次，然後加總抽樣結果，並且畫成圖 8.4 如下：
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 8.4: 丟兩枚硬幣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 8.5: 擲兩枚硬幣1000次的正面次數長條圖二
Figure 8.6: 打珠台
Figure 8.7: 打珠台事件與獎品
8.2.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})\)，就能計算變異數，以表8.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)等等。
8.2.3 二元（白努利）分佈(Bernoulli distribution)
- 白努利分佈假設變數的值為1或是0。假設\(P(X=1)=p\)，那麼\(P(X=0)=1-p\)，PMF寫成：
\[\begin{align*} p(X=x)=p^\cdot (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^1\cdot (1-p)^{1-1}=p\cdot 1=p\)。當\(x=0\)，\(p(X=0)=1\cdot (1-p)=1-p\)。
- 變異數為\(p\times(1-p)\)。因為\(0\le p\le 1\)，所以變異數大於0。
☛已知一疊撲克牌65張中，有35張是鑽石與紅心，請問抽1張撲克牌，抽到紅色的機率多高？
- 已知\(p=\frac{35}{65}\approx0.543\)，代入公式：
\[0.543^1\times(1-0.543)^{1-1} = 0.543\]
- 模擬抽1000次，得到圖8.8：
Figure 8.8: 從65張撲克牌抽到1張紅色牌的機率
- 求出二元分布的期望值：
\[\begin{align} \mathbf{E}[X]& = \sum x\cdot P(X=x)=0\cdot (1-p)+1\cdot p\\ & = p \end{align}\]
- 求出二元分布的變異數：
\[\begin{align} \text{Var}[X]& = \mathbf{E}[X^2]-(\mathbf{E}[X]^2)\\ & = 0^2\cdot (1-p)+1^2\cdot p - p^2 \\ & = p(1-p) \end{align}\]
8.2.4 二項分佈(Binomial distribution)
- 假如已知某個事件為真的機率為\(p\)，我們進行\(n\)次抽樣，其中\(X=x\)次得到某事件為真，該隨機變數服從二項分布。二項分布的機率質量函數寫成：
\[\begin{align} P(X=x)= \binom{n}{x}p^x(1-p)^{n-x}=\frac{n!}{x!(n-x)!}p^x(1-p)^{n-x} \tag{8.2} \end{align}\]
- 要注意的是如果進行\(N\)次實驗，每一次實驗有\(n\)次抽樣，其中\(x\)次是實驗成功的事件，則公式(8.2)\(\times N\)會得到期望的成功事件數。
☛ 假設白天酒駕的機率為\(p=10\%\)，如果警察隨機攔檢20輛，且假設攔檢彼此互相獨立，剛好出現4位酒駕的機率是多少？
答：20次當中出現4次酒駕的機率質量函數為\[P(X=4)=\binom{20}{4}0.1^{4}(1-0.1)^{20-4}\]
可用R的
dbinom(4, 20, 0.1)
## [1] 0.08978
Figure 8.9: 二項分配之次數與機率圖
- 也可以直接用指令計算，注意\(\frac{n!}{(n-k)!}\)的指令為
choose(n,k) ：
p=0.1
n=20
x=4
choose(n, x)*(p^x)*(1-p)^(n-x)[1] 0.08978
#or
factorial(n)/(factorial(x)*factorial(n-x))*(p^x)*(1-p)^(n-x)[1] 0.08978
- 二項分布的期望值剛好是二元分佈的期望值的\(n\)倍，\(n\)是實驗的次數，所以寫成\(E(X)=np\)。而變異數\(V(X)=np(1-p)\)。
8.2.5 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\)是在一定期間內的事件發生次數的平均數。圖8.10顯示三種平均數的Poisson分佈：
Figure 8.10: 3種平均數的Poisson分佈
可以發現，當Poisson分佈的平均數\(\lambda\)越大，越接近常態分佈。
如果進行n次白努利試驗，結果會近似Poisson分佈。這是因為二項分佈的期望值是\(n\times p\)，剛好是Poisson分佈的期望值\(\lambda\)。
當\(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.047168.2.6 Poisson分佈實例1
現代統計學 第197頁，已知平均5分鐘有1.5人使用ATM，半小時平均有9人，請問半小時有5人使用ATM的機率是多少？- 我們可以應用
dpois(k, \(\lambda\)) ，求出當平均值\(\lambda\)為9的時候，x=5的機率：
lambda=9
dpois(5, lambda)[1] 0.06073
半小時內有6個人來用ATM的機率只有6%。
或者用公式計算：
lambda=9
x = 6
p=lambda^{x}*exp(-lambda)/factorial(x)
cat("P(X=6)=", p)P(X=6)= 0.09109 - 圖??表示\(\lambda=9\)時的PMF：
# Set lambda
lambda <- 9
# Define the range of x values (usually around lambda ± 3*sqrt(lambda))
x_vals <- 0:20
# Calculate PMF values
pmf_vals <- dpois(x_vals, lambda = lambda)
# Plot
barplot(pmf_vals,
names.arg = x_vals,
col = "skyblue",
main = expression(paste("Poisson PMF with ", lambda, " = 9")),
xlab = "Number of Events (x)",
ylab = "Probability",
border = "blue")
- Poisson 分佈的累積機率密度函數如圖8.11：
# Define average rate (lambda)
lambda <- 6 # 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 = 9") +
theme_bw()
p + ggplot2::geom_vline(xintercept = 6, linetype="dashed", color = "#DD1100")
Figure 8.11: Poisson 分佈的累積機率密度函數
- Poisson分佈的PDF是\(P(X=k, \lambda)\)，把k對應的機率累積，就可以得到圖 8.11 顯示的Poisson分佈的累積機率。\(f(X)=P(X>k|\lambda)\)。當\(k\)越大，\(P(X>k|\lambda)\)也越大，然後就接近1。曲線上的每一個點代表當X=k時，函數所累積的機率。
8.2.7 Poisson分佈實例2
現代統計學 第202頁，已知咖啡店平均每小時有50位顧客，請問每小時需要幾位服務人員，才不會讓每位服務人員需要服務超過4位客人的機率高於0.5？換句話說，\(\lambda\)應該等於多少，才不會出現一個服務生必須服務超過4個客人？也就是\(k>4\)。如果用\(b\)表示服務生人數，\(\lambda=50/b\)。我們要求\(P(X\le k\mid \lambda)\)。
圖 8.12顯示機率密度函數， 當平均服務人數\(k\)從0增加到20時，\(P(X=k\mid \lambda =5)\)在\(k = 5\)時達到最高點，也就是\(P(X=5\mid \lambda=5)\)接近0.175，而\(P(X=10\mid \lambda=5)\approx 0.018\)。當有10位服務人員時，按照過去經驗，每小時來50人，所以每一位很可能需要服務5人，機率是0.018。但是我們有興趣的是每一位工作人員必須服務超過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 8.12: Poisson機率分佈之機率密度圖
- 圖 8.13顯示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\)時，也就是工作人員10人，平均服務人數超過4人的機率仍然有55.9%，這是因為過去平均每小時有50位顧客。
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 8.13: 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計算從曲線另一邊累積的機率。圖 8.14 顯示當b > 10，\(P(X\le 4\mid 5)\)的累積機率低於0.5，因此，咖啡廳應該雇用至少11位服務人員，才會讓一位服務人員需要服務超過4位客人的機率低於0.5。
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 8.14: 超出特定平均值之累積機率密度圖
9 連續變數
- 連續變數的定義為：如果有一個函數，對任何一個區間，也就是\(A\subset R\)，有如下的關係就是連續變數：
\[ P(X\subset A)=\int _{A}f_{x}(x)dx \]
其中，\(f_{x}(x)\)稱為機率密度函數，也就是Probability Density Function(PDF)。
連續變數的公設有：
- \(f(x)\ge 0\quad \forall\hspace{.2em}x\) (non-negativity)
- \(\int_{-\infty}^{\infty} f(x)dx=1\) (Total probability =1)
- \(P(A)=P(a\le X\le b)=\int_{A} f_{x}(x)dx\) (Additivity)
9.1 常態分佈的公設
- 圖9.1顯示\(f(x)>0\)，而且曲線下面積\(\approx 1\)，陰影的機率是\(P(-1\le X\le 1)\)-1到1的積分\(\int_{-1}^1 f(x)dx\)。
Figure 9.1: 常態分佈
9.2 連續變數與機率分佈
連續變數常用的機率分佈有：常態分佈、標準常態分佈、均等分佈、指數分佈等等。
定義：連續隨機變數\(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}\text{exp}\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\]
- 隨機變數X如果是常態分佈，寫成：
\[X\sim N(\mu, \sigma^2)\]
- 我們可以透過\(Z=\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}\]
9.2.1 連續與不連續變數圖形：相對次數與相對次數密度
- 例如，我們從常態分佈中抽出1,000個值，然後畫兩個長條圖如圖 9.2。第一個圖是把1,000個值分成若干組，圖一顯示相對次數;圖二則是相對次數密度。相對次數=相對次數密度乘以組距。
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="相對次數密度")
Figure 9.2: 單位與實數組距的長條圖
- 表9.1可以幫助理解：
library(kableExtra)
d<-rbind(head=c('Class Interval','Frequency','Relative Frequency', 'Class Width', 'Relative Frequency Density'),
X=c('0-10','20', '0.20','10','0.20/10=0.02'),
Y=c('10-30', '30', '0.30', '20', '0.30/20=0.015'),
Z=c('30=50', '50', '0.50', '20', '0.50/20=0.025'))
row.names(d) <- c('','','','')
kbl(d,
booktabs=T, caption = '100個觀察值分成三組') %>%
kable_paper("striped", full_width = F) %>%
row_spec(1, bold = T, background = "gray65")| Class Interval | Frequency | Relative Frequency | Class Width | Relative Frequency Density | |
| 0-10 | 20 | 0.20 | 10 | 0.20/10=0.02 | |
| 10-30 | 30 | 0.30 | 20 | 0.30/20=0.015 | |
| 30=50 | 50 | 0.50 | 20 | 0.50/20=0.025 |
9.2.2 機率與機率密度函數
- 假設我們相信\(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\]
- 例如常態分佈的部分面積如圖 9.3：
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 9.3: 機率分佈曲線下的區域
- 我們可以用積分求0到1之間曲線底下的面積，得到\(0\le x\le 1\)的機率，也就是套用積分在常態分佈的函數：
\[f(x)=\frac{1}{\sqrt{2\pi}\sigma}\text{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個標準差（或者Z值），那麼面積就是\(68\%\)。
累加機率函數(cumulative density function, CDF)代表機率從F(X=\(-\infty\))一直累加到F(X=x)，表示為：
\[F(X)=P(X\leq x)=\int_-\infty^xf(t)dt\]
- 以圖9.4表示常態分佈而且\(P(X\leq 1)\)的CDF：
# Define the range for the plot
x <- seq(-3, 3, length = 100) # Range of x from -3 to 3
# Compute CDF values for a standard normal distribution
cdf_values <- pnorm(x, mean = 0, sd = 1)
# Plot the CDF curve
plot(x, cdf_values, type = "l", col = "blue", lwd = 2,
main = "CDF of Normal Distribution with Shaded Area",
xlab = "x", ylab = "P(X ≤ x)", ylim = c(0, 1))
# Shade the area under the CDF curve from -∞ to x = 1
#polygon(c(-3, x[x <= 1], 1), c(0, cdf_values[x <= 1], 0),
# col = "#ffff11", border = NA)
# Add a vertical line at x = 1
abline(v = 1, col = "red", lty = 2)
# Add a horizontal line to intersect x = 1
abline(h = 1 - pnorm(0, 1), col = "red", lty = 2)
text(x = -2, y = 0.75, labels = paste0("P(X ≤ 1) ≈ ", round(1 - pnorm(0, 1), 4)), cex = 1.2, col = "black")
# Add grid
grid()
Figure 9.4: 常態分佈的累積機率密度
- 以圖9.5表示常態分佈而且\(P(X\leq 1)\)的PDF：
# Define the range of x values for the plot
x <- seq(-3, 3, length = 1000) # Range of x from -3 to 3
# Compute PDF values for the standard normal distribution
pdf_values <- dnorm(x, mean = 0, sd = 1)
# Compute the cumulative probability (shaded area) P(X ≤ 1)
area_prob <- pnorm(1, mean = 0, sd = 1)
# Plot the PDF curve
plot(x, pdf_values, type = "l", col = "blue", lwd = 2,
main = "PDF of Standard Normal Distribution with Shaded Area",
xlab = "x", ylab = "Density", ylim = c(0, 0.5))
# Shade the area under the curve from -∞ to x = 1
x_shade <- seq(-3, 1, length = 500) # Define x values for shading
y_shade <- dnorm(x_shade, mean = 0, sd = 1) # Corresponding y values
polygon(c(-3, x_shade, 1), c(0, y_shade, 0), col = "#bbee11", border = NA)
# Add a vertical line at x = 1
abline(v = 1, col = "red", lty = 2, lwd = 2)
# Label the shaded area probability
text(x = -0.4, y = 0.1, labels = paste0("P(X ≤ 1) ≈ ", round(area_prob, 4)), cex = 1.2, col = "black")
Figure 9.5: 標準常態分佈的PDF以及機率
假設c<d，則F(c\(\le\)X\(\le\)d)=F(d)-F(c)
對cdf進行微分可得機率密度函數表示為：
\[F^{`}(x)=f(x)\]
\[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\\ F_{x}(x)=\frac{dF(x)}{dx}=f_{x}(x)^{-\infty} \]
- a到b之間的CDF的期望值為：
\[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}\]
9.2.3 不連續函數的累積機率密度應用
- 參考林惠玲與陳正倉(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}\]
9.3 期望值與變異數
9.3.1 期望值
\[\boxed{E(X)=\int_{a}^{b}xf(x)dx=\mu} \]
- 在間斷變數的例子，變數有個別的變量，例如變數\(X\)有四個可能的值1, 2, 3, 4，發生的機率如表 9.2：
| 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 9.6: 間斷變數之相對次數圖
\[\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*}\]
- 圖 9.7 顯示四種大小不同的\(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 9.7: 連續變數之機率密度圖
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 9.8: 機率密度圖的區域
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計算。
9.3.2 變異數
- 在了解連續變數的平均數之後，我們想要知道變數的離散程度。連續變數的變異數的計算方式為：
\[\begin{align} V(X)=\int_{a}^{b}(x-\mu)^2\cdot f(x)dx \tag{9.1} \end{align}\]
- 方程式 (9.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{9.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)\)的期望值以及變異數？
9.4 常態機率分佈
- 常態分佈的值落在平均數加減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 9.9: 常態機率分佈
- 如果資料近似鐘形的常態分佈，可以在長條圖上面加上常態分佈曲線，如圖 9.10：
#抽樣然後加總得到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 9.10: 近似常態分佈之長條圖加上常態分佈曲線
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 9.11: 累積機率密度圖
R提供的函數畫累積密度曲線，如圖 9.12。當\(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 9.12: R函數之累積密度曲線圖
- 圖 9.13 顯示當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 9.13: 機率分佈曲線下的一半區域
9.4.1 常態分配相關函式
R的常態分佈函式有dnorm 產生機率密度函數，也就是代入X之後的值。我們以IQ這個連續變數為例，最小為70、最高為150，平均值110、標準差則是任意的15。我們可以產生對應的機率密度，並且畫成圖9.14：
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 9.14: IQ機率密度圖
pnorm 產生累積機率密度函數，也就是代入X之後產生從最小值一直累積到X的機率。畫成圖 9.15：
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 9.15: IQ累積機率密度圖
pnorm 產生累積機率密度函數，qnorm 則是pnorm 的相反，也就是產生從最小值一直累積到X的機率所對應的值。畫成圖9.16：
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 9.16: 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 產生來自常態分佈的隨機亂數，我們產生三個不同樣本規模但是相同的平均數以及變異數的樣本，然後畫相對次數的長條圖 9.17，可以發現，樣本規模越大，則越集中在平均數附近。
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 9.17: 三種樣本規模的長條圖
9.5 均等機率分佈
- 均等分配指隨機變數在某一連續的區間，有同等的機率密度，例如等車時間、飛行時間、排隊買美食時間。可以寫成:
\[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{9.3} \end{equation*}\]
- 在方程式(9.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)，\)圖 9.18 顯示均等機率分佈曲線下的區域：
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 9.18: 均等機率分佈曲線下的區域
- 由以下可知，\(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
- 圖 9.19 顯示日出時間是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 9.19: 一致機率分佈曲線下的區域
9.5.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{9.4} \end{align*}\]
\[\begin{equation*} f(x) = \begin{cases} \frac{1}{X} & 0\le x\le 1\\ 0 & \text{otherwise} \end{cases}\tag{9.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\]
- 或者根據方程式(9.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
9.5.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 9.20: 均勻分佈的累積機率圖
qunif 回傳均勻分佈的百分位如圖9.21：
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 9.21: 均勻分佈的百分位圖
runif 產生均等分配的隨機變數，例如抽出位於(0,10)的1000個數字，繪成長條圖9.22之後，每一個數字的機率密度約為0.1：
set.seed(1000)
rand.unif<-runif(1000, 0, 10)
hist(rand.unif, freq = FALSE, col='#4433ae', density = 40, xlab='x', main='')
Figure 9.22: 均等分配隨機變數的長條圖
punif(1, 0, 10)[1] 0.1
9.6 指數機率分佈
\[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}\)。
9.6.1 指數分配期望值與變異數
\[E(X)=\frac{1}{\lambda}\]
\[V(X)=\frac{1}{\lambda^2}\]
9.6.2 指數分配的相關指令
dexp 回傳機率密度，以上一題為例，可假設\(0\leq X\leq 50\)，得到對應的機率密度由高而低如圖9.23：
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 9.23: 指數分配機率密度圖
- 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\)，得到對應的機率密度由高而低如圖9.24：
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 9.24: 指數分配累積機率密度圖
- 圖9.25顯示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 9.25: 指數分配累積機率區域圖
- 我們可以用
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}\)由低而高如圖9.26：
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 9.26: 指數分配百分位圖
- 最後，
rexp 回傳指數分配\(\lambda\)值的隨機亂數，可排列成為機率密度圖9.27：
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 9.27: 指數分配機率密度圖
- 如果計算\(\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張訂單。
10 機率作業
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，這兩個事件是否為獨立事件？11 隨機變數作業
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}\]
12 更新講義時間
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)最後更新時間: 2025-04-16 20:27:13