setwd("D:/")
school <- read.csv("台灣地區大專院校學生統計表.csv")

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
school %>% filter(等級別 %in% c("B 學士","M 碩士","B 四技")) -> new
q <- NULL
for(i in 1:500){
  n = 5
  a=sample(1:148,n,replace = F)
  aa=unique(new$學校名稱)[a]
  
  k = new %>% filter(學校名稱%in%aa) %>% group_by(學校名稱)  %>% mutate(sum = sum(總計))
  k$sum %>% unique() -> a 
  sp2 <- (sum(a^2)/4)/(n-1)
  D <- (0.15^2)*(6292.6^2)/4
  
  q[[i]] <- 148*sp2/((148*D)+sp2)
}
hist(q , breaks = 10)