以學測篩選倍率為例-介紹統計模擬並用R語言實作

詹永裕
詹永裕
Nov 8 · 15 min read

沒有真實的資料還是能做分析! -統計模擬簡介

做資料分析最大的問題就是...沒有資料! 沒有資料要怎麼做分析?此篇文章介紹統計模擬的概念並以學測篩選倍率為例,讓你即使沒有真實資料也能夠做分析!

相信各位都有經歷過高三那段日夜苦讀拼學測、考大學的時光,那自然會明白何謂篩選倍率以及其重要性,我們知道篩選倍率非常重要,影響了你能不能夠被選進該科系也影響了各科系所篩選進的學生成績分佈,然而篩選倍率實際上是如何影響卻鮮少有文章著墨、探討,因此在這邊介紹之前受成大統計的張升懋教授指導使用統計模擬的方法並用R語言實作探討不同篩選倍率實際上所造成的影響。

http://45.33.36.58/article/detail?id=1087

以學測篩選倍率為例-介紹統計模擬並用R語言實作預計會分成多篇文章,這一篇先跟各位介紹基本的觀念,及如何在R實作,接下來還會比較使用R或Julia執行速度上的差異。

可能有很多人不清楚什麼是統計模擬,那究竟什麼是統計模擬?

先附上一段定義,以下擷取自MBA智庫

蒙特卡羅方法又稱統計模擬法、隨機抽樣技術,是一種隨機模擬方法,以概率和統計理論方法為基礎的一種計算方法,是使用隨機數來解決很多計算問題的方法。將所求解的問題同一定的概率模型相聯繫,用電子電腦實現統計模擬或抽樣,以獲得問題的近似解。

簡單來說,統計模擬就是使用電腦產出隨機數以幫助我們解決一些複雜的計算問題及在沒有實際資料的狀況下模擬實際資料。在此由於我們無法取得學測篩選的真實資料,因此我們使用統計模擬去模擬實際資料。

期望這篇文章能夠讓各位在做科系選擇時更加清楚篩選倍率的影響,同時也幫助各科系設定篩選倍率,更讓對於統計模擬&R有興趣之同學一起學習。

若有對於學測篩選的機制如何運作還不清楚的同學可以先參考這篇文章

以下正文開始:

1. 資料生成

首先設定我們希望產生多少筆資料去做模擬,在此設定為10萬筆。

接著考量到每一個學生(每一筆資料)不同科目的成績彼此會有所影響,例如:若學生若英文好那他在同樣屬於文科的國文表現可能也不錯,反之一個學生若在數學表現較佳那他在同屬於理科的自然有較高的可能也表現優異,因此我們設定一個共變異數矩陣以模擬各科系之間的影響幫助我們生成資料。

# number of students
n = 100000
# C: Chinese; E: English; M: Math; N: Nature; S: SocialC_E = .6; C_M = -.2; C_N = -.2 ;C_S = .4
E_M = -.4; E_N = -.2; E_S = 0
M_N = .8; M_S = 0
N_S = -.2
# 共變異數矩陣
Cor <- matrix(c(1, C_E, C_M, C_N, C_S,
C_E, 1, E_M, E_N, E_S,
C_M, E_M, 1, M_N, M_S,
C_N, E_N, M_N, 1, N_S,
C_S, E_S, M_S, N_S, 1), 5, 5)

接下來定義兩個函數,分別是DGP以及Cal_Grade_Level,DGP幫助我們用常態分配和共變異數矩陣生成每個學生各科的原始分數,再使用Cal_Grade_Level這個函數幫助我們把原始分數轉成大家所熟知的級分,為了方便,我們直接把Cal_Grade_Level包進DGP了。

在這邊最重要的概念是使用常態分配和共變異數矩陣生成資料了,雖然真實資料並不一定是常態分配,但我們無法取得真實資料因此假設資料為常態去模擬以利接下來的分析,這就可以算是統計模擬,生成隨機數來模擬真實的狀況,最重要的是這樣的假設也不會太不合理以至於影響後續結果。

head(DGP(n, Cor))Chinese English Math Nature Social
[1,] 9 12 5 7 9
[2,] 5 7 9 9 7
[3,] 8 11 9 9 5
[4,] 9 11 5 7 7
[5,] 7 7 7 9 8
[6,] 11 12 9 12 8

利用以上的函數在將人數、共變異數矩陣放入便可產出以下的資料,每一列代表一位學生,每一個欄位代表各學生該科的成績。

DGP以及Cal_Grade_Level定義如下

# 學測及分算法可參考這篇文章Cal_Grade_Level <- function(x){
h <- quantile(x, .99)
h <- round(mean(x[which(x >= h)])/15, digit=2)
cutpoint <- h*(1:15)
grade <- rep(0, length(x))
for(i in 1:15) grade <- grade + ifelse(x > cutpoint[i], 1, 0)

return(grade)
}
DGP <- function(n, Cor){X <- matrix(rnorm(n*5), n, 5)%*%chol(Cor)
# 用Cor 和 rnorm 生成 n X 5的矩陣
X <- X - min(X); Z <- X
for(i in 1:5) Z[,i] <- Cal_Grade_Level(X[,i]) # 轉換成級分
colnames(Z) <- c("Chinese", "English", "Math", "Nature", "Social")
return(Z)
}

現在我們就有辦法生成學生的成績了~

也就是說已經完成統計模擬中生成資料的步驟了,接下來要做的就是要看看如何分析這筆隨機產出的資料!

2. 篩選檢定標準

生成完資料後就可以進入篩選的程序了~

第一步我們需要先篩選檢定標準,也就是各科系所對各科設定的最低標準「頂標」、「前標」、「均標」、「後標」、「底標」。

所以在這邊我們同樣定義兩個函數Cal_Standard以及Screening_by_STD,Cal_Standard做的是計算生成資料的五標,而Screening_by_STD做的就是給定該科系的篩選標準,回傳通過篩選的學生成績,另外為了方便Cal_Standard直接被我們包進Screening_by_STD了

Score <- DGP(n, Cor) #隨機產出學生原始成績並將原始成績轉為級分STD_vec <- c(3,2,2,0,0) # 給定檢定標準 
# 數字代表1:頂標 2:前標 3:均標 4:後標 5:底標
# 而依順序為 國文、英文、數學、自然、社會
# 因此我們給定的檢定標準為 國文均標、英文前標、數學前標
X_sc <- Screening_by_STD(Score, STD_vec) # 回傳通過給定篩選標準的資料head(X_sc)Chinese English Math Nature Social
[1,] 10 11 12 13 9
[2,] 11 10 11 10 11
[3,] 9 10 10 11 9
[4,] 12 11 11 10 11
[5,] 12 11 10 8 11
[6,] 9 11 10 10 9
nrow(X_sc) # 有多少學生通過篩選?[1] 8227

Cal_Standard以及Screening_by_STD定義如下

Cal_Standard <- function(X){
# input: grade 0~15
out <- sapply(1:5, function(i) quantile(X[,i], c(.88, .75, .50, .25, .12)))
colnames(out) <- colnames(X)
return(out)
}
Screening_by_STD <- function(Score, STD.vec){
STD <- Cal_Standard(Score) # 計算五標
STD_vec <- rep(NA,length(STD.vec))
o <- 1
for (h in 1:length(STD.vec)){
if (STD.vec[h] != 0){
STD_vec[o] <- c(STD[STD.vec[h],STD.vec[o]])}
else{
STD_vec[o] <- 0}
o <- o+1
}

G <- matrix(1, nrow(Score), ncol(Score))
for(i in 1:5) if(STD_vec[i] > 0) G[,i] <- ifelse(Score[,i] >= STD_vec[i], 1, 0)
G <- apply(G, 1, sum)
SC_X <- Score[which(G == 5),]
colnames(SC_X) <- colnames(Score)
return(SC_X)}

以上我們就完成了生成資料以及檢定標準的篩選了,我們對隨機生成的資料進行檢定標準的篩選,接下來就可以進入我們最重要的部份~ 篩選倍率!

3. 篩選倍率

完成了篩選檢定標準,就到了篩選倍率的階段了。

在這邊我們同樣定義一個函數Selection幫助我們進行篩選倍率,Selection所做的是給定篩選倍率和篩選人數,回傳通過篩選的學生。

以下的程式碼可以看到,給定篩選倍率 有哪些學生通過篩選倍率能進入第二階段面試。

SEL <- rbind(c(0,1,0,0,0,9),
c(0,0,1,0,0,3)) # 給定篩選倍率
# 同樣的 vector依序為 國文、英文、數學、自然、社會 不同的是最後面代表的是篩選的倍率
# 因此我們給定的篩選倍率為 英文9倍 數學3倍
out <- Selection(X_sc, SEL, 14)head(out)
Chinese English Math Nature Social
[1,] 14 15 13 13 9
[2,] 14 15 13 12 12
[3,] 13 14 13 14 9
[4,] 14 14 13 13 11
[5,] 12 15 12 13 11
[6,] 12 15 12 11 9

走到這裡我們就完成大部份的模擬階段啦!我們可以使用以下程式碼去看篩選近來的學生之成績分佈。

round(apply(out, 2, mean), digits=2) # 學生五科平均成績Chinese English    Math  Nature  Social 
12.45 14.26 11.36 11.93 9.05
round(apply(out, 2, sd), digits=2) # 學生五科成績之標準差Chinese English Math Nature Social
1.63 0.45 0.76 1.24 2.16

Selection定義如下

Selection <- function(X, SEL, number){

SEL <- matrix(SEL,ncol=6)
n <- nrow(SEL)
if (n != 1){
SEL <- SEL[order(SEL[,ncol(SEL)],decreasing=TRUE),]}
Z <- X
for(i in 1:n){
nn <- ceiling(14*SEL[i,6])
temp <- Z%*%SEL[i,-6]
Z <- Z[order(temp,decreasing = TRUE),]
Z <- Z[1:nn,]
}
return(Z)}

4. 小結

經過以上三個步驟 資料生成、檢定標準、篩選倍率

我們可以將模擬的程式碼整合如下,如此一來我們變能夠使用統計模擬生成給定檢定標準及篩選倍率且最後能夠通過層層篩選的學生之成績了!

# number of students
n = 100000
# C: Chinese; E: English; M: Math; N: Nature; S: SocialC_E = .6; C_M = -.2; C_N = -.2 ;C_S = .4
E_M = -.4; E_N = -.2; E_S = 0
M_N = .8; M_S = 0
N_S = -.2
Cor <- matrix(c(1, C_E, C_M, C_N, C_S,
C_E, 1, E_M, E_N, E_S,
C_M, E_M, 1, M_N, M_S,
C_N, E_N, M_N, 1, N_S,
C_S, E_S, M_S, N_S, 1), 5, 5)
Score <- DGP(n, Cor) # 生成資料
STD_vec <- c(3,2,2,0,0) # 給定檢定標準
X_sc <- Screening_by_STD(Score, STD_vec) # 用檢定標準進行第一次篩選
SEL <- rbind(c(0,1,0,0,0,9),
c(0,0,1,0,0,3)) # 給定篩選倍率
out <- Selection(X_sc, SEL, 14) # 回傳通過篩選倍率的42(14*3)位學生之成績
round(apply(out, 2, mean), digits=2) # 看篩選近來之學生的成績分佈
round(apply(out, 2, sd), digits=2)

5. 還沒結束…

做到這裡當然還沒結束,你還有發現什麼問題嗎?

沒錯!就是 經過我們這樣子的模擬,我們是假設隨機生成的每一位學生都會報名該科系,因為不論是檢定標準、篩選倍率我們都是直接把全部學生的資料都放進去篩選,但這在實際上並不可能!為了更貼近現實的狀況,我們在這邊又多加了一些修改。

首先,我們先將前面的學生剃除,在這邊定義了一個dele的函數,dele幫助我們將總分前百分之幾percent的學生刪掉,因為這些學生可能有更多好的選擇。

Score <- DGP(n, Cor)
round(apply(Score, 2, mean), digits=2)
Chinese English Math Nature Social
9.13 9.32 9.31 9.31 9.13
Score <- dele(Score,0.1) # 刪去前百分之10的學生
round(apply(Score, 2, mean), digits=2)
Chinese English Math Nature Social
8.88 9.16 9.13 9.12 8.95

經過dele的函數可以發現,前百分之10的學生被排除在外了,因此整體學生的平均成績下降。

dele函數定義如下

dele <- function(Score,percent){
n <- nrow(Score)
Score <- Score[order(rowSums(Score),decreasing = TRUE),]
Score <- Score[(n*percent):n,]
return(Score)
}

除了這樣還有一件事沒被考量到,那就是並不是所有學生都會選填該科系,但我們也無法得知實際上到底有哪些學生會選填該科系,因此我們使用隨機抽樣,在有能力經過檢定標準的學生裡選出n個去選填該科系。

在這邊我們定義了samp函數,在所有學生裡選出n位學生申請該科系。

Score <- DGP(n, Cor)
Score <- dele(Score,0.1) # 刪去前百分之10的學生
STD_vec <- c(3,2,2,0,0) # 給定檢定標準
SEL <- rbind(c(0,1,0,0,0,3),c(0,0,1,0,0,6))
# 給定篩選倍率
X_sc <- Screening_by_STD(Score, STD_vec) # 篩選通過檢定標準的學生res <- samp(X_sc,SEL,100,14)
# 在所有通過檢定標準的學生裡隨機抽100為學生選填該科系

samp做的就是取代了Selection函數的功能,而多增加了要抽多少人去選填該科系。

samp定義如下

samp <- function(X_sc,SEL,x,n){
res <- matrix(rep(NA,5), ncol = 5)
X_sc <- X_sc[sample(nrow(X_sc), x), ]
res <- Selection(X_sc, SEL, n)
return(res)
}

如此一來,我們就有了刪去前面幾%的學生,也使用隨機抽樣保證不是所有學生都會選填該科系的問題。

6. 總結

快結束了!能夠走到這邊很不容易,以上重重的步驟可以寫成以下

# number of students
n = 100000
# C: Chinese; E: English; M: Math; N: Nature; S: SocailC_E = .6; C_M = -.2; C_N = -.2 ;C_S = .4
E_M = -.4; E_N = -.2; E_S = 0
M_N = .8; M_S = 0
N_S = -.2
Score <- DGP(n, Cor) # 生成資料STD_vec <- c(3,2,2,0,0) # 給定檢定標準
X_sc <- Screening_by_STD(Score, STD_vec) # 用檢定標準進行第一次篩選
SEL <- rbind(c(0,1,0,0,0,9),
c(0,0,1,0,0,3)) # 給定篩選倍率
res <- samp(X_sc,SEL,100,14) # 隨機抽選100位學生申請該科系並回傳通過篩選倍率的42(14*3)位學生之成績round(apply(res, 2, mean), digits=2) # 看篩選近來之學生的成績分佈
round(apply(res, 2, sd), digits=2)

如此一來我們就能夠隨機生成資料和給定檢定標準、篩選倍率以及依照該科系去客製化要刪去前幾百分比的學生和會有多少人申請的學生成績分佈。這樣我們便可以使用各式各樣不同組合的篩選倍率去觀察篩選倍率實際上會有什麼影響。

有興趣的同學也可以自行使用這些程式碼去模擬感興趣的科系所設定的篩選倍率會篩選進怎麼樣的學生,怎麼樣的學生才能被篩選進來?

學測篩選倍率探討就是一個統計模擬的例子,雖然我們沒有真實的資料,但我們依然能夠利用一些技巧去隨機生成資料並加以分析嘗試去模擬真實的狀況。

今天的介紹到這邊就告一段落,謝謝!

題外話 宣傳一下 Julia's Knights

https://www.facebook.com/juliasknights/

這是我們成大統計系學生一起發起的活動,希望大家多支持,按讚我們的粉專,喜歡我的文章的話也可以幫我拍手: )

詹永裕

Written by

詹永裕

統計、資料分析、分享生活

Welcome to a place where words matter. On Medium, smart voices and original ideas take center stage - with no ads in sight. Watch
Follow all the topics you care about, and we’ll deliver the best stories for you to your homepage and inbox. Explore
Get unlimited access to the best stories on Medium — and support writers while you’re at it. Just $5/month. Upgrade