2013.12 貝氏應用統計期末報告，附錄程式碼

0_Readme.md

Readme

這些程式碼包含：

• 資料前處理
• 繪圖
• Bayesian linear regression model 計算

三大部份的程式碼，以下將分別列述。 由於原始資料需經過一定處理轉換才能使用，但細節較多且複雜，另獨立列述。

如果要跳過資料前處理，可以使用使用 101_境外學生總計_國別_ConsumerIndex.csv link 作後續的繪圖與 model 計算。

所有的檔案可以一次下載，使用左方的 Download Gist 一連結。

原始資料來源

• 101_大專校院境外學生在臺人數統計，教育部統計處 http://stats.moe.gov.tw/files/detail/101/101_ab101.xls 為求清楚、方便處理，重新命名並用新的 Excel 格式另存為 101_大專校院境外學生在臺（國別、校別）.xlsx

• Cost of Livings Index by Country, Numbeo 這邊使用 HTML 直接複製原始 102 個國家表格於 Excel 中，另存為 CSV 格式 living_cost.csv，見本 gist 後面檔案。 但由於許多國家的資料在 Numbeo 分別瀏覽時有，但在此表中沒有出現，故另行補上

以上的資料在處理時放在 dataset/ 資料夾中，程式碼的相對路徑皆需要根據使用者設定作修改。

資料前處理

包含兩檔案：

• parse_raw_data.ipynb (open with nbviewer)
• living_cost.ipynb (open with nbviewer)

使用 IPython Notebook 撰寫，檔名後連結為一個線上觀看 notebook 格式的網站。 裡面所儲存的程式碼使用 Python 語言撰寫。如果需要在本機端開啟並重新執行，需要安裝 Python 3.3.2 以及相關的套件：

• Pandas (提供向 R 的 data.frame 功能）
• Numpy 1.7+
• IPython 1.0+
• pygeocoder
• geopy

安裝方式

在非 Windows 系統上安裝，以 Ubuntu 舉例

# Python 3.x dependcies
sudo apt-get install build-essential git xz-utils \
  libncursesw5-dev libreadline-dev libssl-dev libgdbm-dev \
  libc6-dev libsqlite3-dev tk-dev

# Install Python 3.3.2
cd /tmp
wget http://python.org/ftp/python/3.3.2/Python-3.3.2.tar.xz
tar Jxvf Python-3.3.2.tar.xz
cd Python-3.3.2
./configure
make
sudo make install

# Install setuptools
wget https://bitbucket.org/pypa/setuptools/raw/bootstrap/ez_setup.py
sudo python3 ez_setup.py
sudo easy_install-3.3 pip

# Install Numpy
sudo apt-get install liblapack-dev libatlas-base-dev gfortran
sudo pip-3.3 install numpy

# Install IPython
sudo easy_install-3.3 ipython[all]

其餘套件均可以使用 sudo pip-3.3 install <package_name> 方式安裝。 就可以使用 ipython notebook 執行相關的程式碼

繪圖

繪圖的程式碼較複雜，故另外獨立出來。主要有兩程式碼

• display_distance.R : 距離相關的繪圖
• display_livingcost.R : 消費指數相關的繪圖

需要安裝像 ggplot2, gridExtra, googleVis 等套件。

Model 計算

• Full model (把 4 個 X 放入計算): full_linear_reg_model.R
• Model selection: model_select.R

101_境外學生總計_國別_ConsumerIndex.csv

display_distance.R

##########################
# 留學生人數與距離的關係 #
##########################

# 資料前處理 ########
df.bycountry <- read.csv("dataset/101_境外學生總計_國別.csv", stringsAsFactors=FALSE)
colnames(df.bycountry) <- c("continent", "country_orig", "student_total", 
                            colnames(df.bycountry)[-(1:3)])
# 去除 大陸、香港、澳門
# df.bycountry <- df.bycountry[-c(2, 5, 8), ]
# 去除 大陸
# df.bycountry <- df.bycountry[-c(8), ]

# 全部國家的繪圖 #######
library(ggplot2)
g.full <- ggplot(df.bycountry, aes(x=tpe_distance, y=student_total, color=continent))
# g + geom_point() + scale_y_log10() + scale_x_log10()  # 使用 log-scale 作圖
g.full <- g.full + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  annotate("text", x=5000, y=8500, label="馬來西亞", family="Heiti TC Medium", size=6) +
  annotate("text", x=12000, y=4500, label="美國",family="Heiti TC Medium", size=6) +
  annotate("text", x=3000, y=17500, label="中國大陸",family="Heiti TC Medium", size=6) +  
  theme(legend.position=c(1,1), legend.justification=c(1,1))

## 取前 50% 多的國家 #######
g.upper_median <- ggplot(df.bycountry[df.bycountry$student_total > 30, ], 
                         aes(x=tpe_distance, y=student_total, color=continent))
g.upper_median <- g.upper_median + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  theme(legend.position=c(1,1), legend.justification=c(1,1))

## 取另外一半 #######
g.lower_median <- ggplot(df.bycountry[df.bycountry$student_total <= 30, ], 
                         aes(x=tpe_distance, y=student_total, color=continent))
g.lower_median <- g.lower_median + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  theme(legend.position=c(1,1), legend.justification=c(1,1))

## 全部合併
library(gridExtra)
grid.arrange(
  arrangeGrob(g.full, ncol=1),
  arrangeGrob(g.upper_median, g.lower_median, ncol=1),
  ncol=2,
  sub=textGrob("與臺灣距離 (km)\n", gp=gpar(fontsize=20, fontfamily="Heiti TC Medium")),
  left=textGrob("\n留學生人數", gp=gpar(fontsize=20, fontfamily="Heiti TC Medium"), rot=90)
)

display_livingcost.R

#####################################
# 留學生人數與 ConsumerIndex 的關係 #
#####################################

# 資料前處理
df.living <- read.csv("dataset/101_境外學生總計_國別_ConsumerIndex.csv")
colnames(df.living) <- c("continent", "country_orig", "student_total", "long_name", "short_name", "longitude", "latitude", colnames(df.living)[-(1:7)])
df.living <- na.omit(df.living) # 去除沒有完整資訊的國家

## 作圖於 google Map 上，這是互動網頁，會需要另存成 png
library(googleVis)
plot(gvisGeoMap(df.living, locationvar="long_name", numvar="ConsumerIndex", 
                options=list(width=1200, height=800, dataMode='regions',
                             colors='[0xFDFF73, 0xFF7573, 0xDC0300]')
))

# 全部國家 #######
g.full <- ggplot(df.living, aes(x=ConsumerIndex, y=student_total, color=continent)) 
g.full <- g.full + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  theme(legend.position=c(1,1), legend.justification=c(1,1)) + 
  geom_vline(aes(xintercept=x, color=continent), data=data.frame(x=39.86, continent="亞洲"), 
             size=1, alpha=0.6, linetype="dashed")
g.full

# 比台灣 CI 高的國家 #######
g.upperCI <- ggplot(df.living[df.living$ConsumerIndex > 40, ], 
                    aes(x=ConsumerIndex, y=student_total, color=continent))
g.upperCI <- g.upperCI + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  theme(legend.position=c(1,1), legend.justification=c(1,1)) + 
  geom_vline(aes(xintercept=x, color=continent), data=data.frame(x=39.86, continent="亞洲"), 
             size=1, alpha=0.6, linetype="dashed") + scale_y_log10()
g.upperCI

# 比台灣 CI 低的國家 #######
g.lowerCI <- ggplot(df.living[df.living$ConsumerIndex <= 40, ], 
                    aes(x=ConsumerIndex, y=student_total, color=continent))
g.lowerCI <- g.lowerCI + geom_point(size=5, alpha=0.6) + theme_bw() +
  theme(text=element_text(family="Heiti TC Medium", size=20)) + 
  labs(x=NULL, y=NULL, color=NULL) +
  theme(legend.position=c(1,1), legend.justification=c(1,1)) + 
  geom_vline(aes(xintercept=x, color=continent), data=data.frame(x=39.86, continent="亞洲"), 
             size=1, alpha=0.6, linetype="dashed") + scale_y_log10()
g.lowerCI

## 合併繪圖
library(gridExtra)
grid.arrange(
  arrangeGrob(g.full, ncol=1),
  arrangeGrob(g.upperCI, g.lowerCI, ncol=1),
  ncol=2,
  sub=textGrob("消費含租屋指數(NYC=100, TW=40)\n", gp=gpar(fontsize=20, fontfamily="Heiti TC Medium")),
  left=textGrob("\n留學生人數", gp=gpar(fontsize=20, fontfamily="Heiti TC Medium"), rot=90)
)

full_linear_reg_model.R

#########################################################
# Bayesian Linear Normal Regression Model using g-prior #
#########################################################

# Data Preprocessing ########
df.living <- read.csv("dataset/101_境外學生總計_國別_ConsumerIndex.csv", stringsAsFactors=FALSE)
colnames(df.living) <- c("continent", "country_orig", "student_total", colnames(df.living)[-(1:3)])
df.needed <- df.living[, c("student_total", "tpe_distance", "ConsumerIndex")]
df.needed <- na.omit(df.needed)  # Trim any with NAs
names(df.needed) <- c("y", "x_dist", "x_ci")
row.names(df.needed) <- NULL


# Construct all Xs ########
## Here we have 
##   x_const=1, x_dist, x_ci, x_interact 
## four regressors. For simplicity, we don't consider further classification
## based on the preliminary view of data.
df.needed$x_const <- 1
df.needed$x_interact <- df.needed$x_dist * df.needed$x_ci
Y <- df.needed$y
X <- as.matrix(df.needed[, -c(1)])


# Sketch ########
panel.cor <- function(x, y, digits=2, prefix="", cex.cor) {
  usr <- par("usr"); on.exit(par(usr))
  par(usr = c(0, 1, 0, 1))
  r = (cor(x, y))
  txt <- format(c(r, 0.123456789), digits=digits)[1]
  txt <- paste(prefix, txt, sep="")
  if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
  text(0.5, 0.5, txt, cex = 5)  # `cex` or `cex * abs(r)`
}
pairs(~ x_dist + x_ci + x_interact, upper.panel=panel.cor, data=df.needed)

# Model Computation ########
## using OLS, x_const will be dropped for singularity
model <- lm(y ~ x_dist + x_ci + x_const + x_interact, data=df.needed)
beta_ols <- model$coefficients[c(2, 3, 1, 5)]  # notice the coeff order should match
dim(beta_ols) <- c(1, 4)
# by definition, simga_ols^2 = sum[(Y - predict_y)^2] / (N-p)
simga_ols <- sqrt(sum((Y - beta_ols %*% t(X))^2) / (length(Y) - dim(X)[2]))

## using g-prior with least prior information
S <- 100000
g <- length(Y); nu0 <- 1; sq2_0 <- simga_ols^2  # prior parameters
n <- length(Y); p <- dim(X)[2]  # some useful info about data
inv_Xsq <- solve(t(X) %*% X)
Hg <- ( g/(g+1) ) * 
  X %*% inv_Xsq %*% t(X)  # beta_ols = inv_Xsq %*% t(X) %*% Y
SSRg <- t(Y) %*% (diag(1, nrow=n) - Hg) %*% Y

sq2_vec <- 1 / rgamma(S, (nu0+n)/2, (nu0*sq2_0+SSRg)/2)
Vb <- g * inv_Xsq / (g+1)
Eb <- Vb %*% t(X) %*% Y

E <- matrix(rnorm(S*p, 0, sqrt(sq2_0)), S, p)
beta <- t(t(E %*% chol(Vb)) + c(Eb))

# Plot ########
par(mfrow=c(2,2))
hist(beta[, "x_const"], main=expression(beta[const]))
hist(beta[, "x_dist"], main=expression(beta[dist]))
hist(beta[, "x_ci"], main=expression(beta[ci]))
hist(beta[, "x_interact"], main=expression(beta[interact]))

living_cost.csv

living_cost.ipynb

model_select.R

# data preprocessing
df.living <- read.csv("dataset/101_境外學生總計_國別_ConsumerIndex.csv", stringsAsFactors=FALSE)
colnames(df.living) <- c("continent", "country_orig", "student_total", colnames(df.living)[-(1:3)])
# df.living <- df.living[-c(8), ]  # 去除 大陸

df.needed <- df.living[, c("student_total", "tpe_distance", "ConsumerIndex")]
df.needed <- na.omit(df.needed)  # Trim any with NAs
names(df.needed) <- c("y", "x_dist", "x_ci")
row.names(df.needed) <- NULL

# Construct all Xs
df.needed$x_const <- 1
df.needed$x_interact <- df.needed$x_dist * df.needed$x_ci
Y <- df.needed$y
X <- as.matrix(df.needed[, -c(1)])

p <- 4
tf <- c(TRUE, FALSE)
models <- expand.grid(replicate(p,tf,simplify=FALSE))
names(models) <- NULL
models <- as.matrix(models)
models <- models[-dim(models)[1],]
models <- models[models[, 3] == TRUE, ]


## model selection
model_select <- function(model) {
  # Ex. model = [T, F, F, T]
  pz <- sum(model)
  Xz <- X[, model]
  g <- n <- length(Y)
  inv_Xsqz <- solve(t(Xz) %*% Xz)
  beta_z_ols <- t(inv_Xsqz %*% t(Xz) %*% Y)
  sq2z_0 <- sum((Y - beta_z_ols %*% t(Xz))^2) / (length(Y) - dim(Xz)[2])
  Hgz <- ( g/(g+1) ) * Xz %*% inv_Xsqz %*% t(Xz)
  SSRgz <- t(Y) %*% (diag(1, nrow=n) - Hgz) %*% Y
  - pz/2 * log(1+n) - log(sq2_0 + SSRgz) * (n+1)/2
}

lml.all=apply(models,1,model_select)
results=cbind(lml.all, models)
order=sort(results[,1],index=TRUE,decreasing=TRUE)
results[order$ix,]

df.results <- data.frame(results[order$ix,])

results[order$ix,1] <- results[order$ix,1] - results[order$ix[1],1]
results[order$ix,1] <- exp(results[order$ix,1])
nconstant <- sum(results[,1])
results[,1] <- results[,1]/nconstant
postprob <- results[order$ix,] 
postprob

df.results$postr_z <- postprob[, 1]
df.results <- df.results[, c(2:5, 1, 6)]
names(df.results) <- c("x_dist", "x_ci", "x_const", "x_interact", "log_ml", "postr_z")
df.results <- df.results[, c(3, 1:2, 4:6)]
df.results[, 5:6] <- round(df.results[, 5:6], 2)
write.csv(df.results, file="model_selection.csv", row.names=FALSE)

parse_raw_data.ipynb