ปัญหา Path ของ Rtools

rtools มันจำเป็นสำหรับ r package ที่มี code จากภาษาอื่น อย่าง c/c++ หรือ fortran ในที่จะต้องมีการ compile ระหว่างการติดตั้ง โดยปกติแล้วมัน จะถูกลงไว้ที่

PATH="C:\Rtools\bin;${PATH}" 

ถ้าลงไว้ที่อื่นและต้องการเรียกใช้ใน R ตัวแปรชื่อ BINPREF จะเป็นตัวแปรที่ R จะเรียกหา path ของ rtools ฉะนั้นถ้าจะเรียกใช้ rtools ที่ลงไว้ที่อื่นก็ต้องเซ็ต path ให้กับตัวแปรนี้ เช่น

BINPREF="X:/R/Rtools-3.5/mingw_$(WIN)/bin/" 

** สังเกตว่า BINPREF ใช้ forward slash นะครับ

ติดตั้ง CmdStan บน Windows 10

ตัวอย่างการติดตั้ง CmdStan (2.17.1) บน Windows 10 ครับ ในที่นี้ ผมใช้กับ Rtools 3.5 ครับ

เริ่มจากการที่เราไป download ตัว zip ของ CmdStan แล้วก็ extract ไว้ที่ไหนสักแห่งตามต้องการครับ

จากนั้นก็เปิด Command Prompt ของ Windows ครับ แล้วใช้คำสั่ง cd เพื่อเปลี่ยนไปยัง path ของ CmdStan ที่ extract ไว้ครับ

แล้วก็ทำการสร้างไฟล์ ชื่อ local (ไม่มีสกุลต่อท้าย) ในโฟลเดอร์ที่ชื่อ make ที่อยู่ path ของ CmdStan โดยในไฟล์ local นี้มีสองรรทัดนี้ครับ

CC=g++

CXX=g++

จากนั้นก็พิมพ์ PATH = c:\rtools\mingw_64\bin;c:\rtools\bin;$PATH เพื่อทำการเซ็ตค่า PATH ให้รู้จักตัว compiler ใน Rtools

แล้วก็ทำการพิมพ์ make build เพื่อทำการ compile ตัว CmdStan ได้เลยครับ

เมื่อเสร็จแล้วจะมีข้อความบอกครับ จากนั้นลองทำสอบกับตัวอย่างที่มาด้วยโดยพิมพ์

make /examples/Bernoulli/Bernoulli.exe

เพื่อสร้างตัวโมเดลจากไฟล์ .stan

จากนั้นก็ลองใช้งานโมเดลกับข้อมูลที่มีมาด้วย โดยพิมพ์

examples\Bernoulli\Bernoulli.exe sample data file=examples/Bernoulli/Bernoulli.data.R

 

เผื่อสงสัยกัน เหตุผลที่ผมใช้ Cmdstan แทนที่จะใช้ rstan เพราะผมใช้งานผ่าน Mathematica ครับ

Stan กับ ปัญหา linear regression

จากที่เคยไปโม้ไว้เยอะว่า Stan มันเจ๋งยังไงให้กับเพื่อนร่วมงานฟัง เมื่อวันศุกร์ที่แล้วผมก็เลยนัดโชว์แบบคร่าวๆพร้อมกับสาธิตการใช้งานกับปัญหาง่ายๆอย่าง linear regression ให้เค้าดูกันครับ ก็เลยคิดว่าจะเอาที่ไปโม้ไว้มาใส่ไว้ที่นี่ด้วยเผื่อว่าใครสนใจใช้งาน Stan แต่ไม่รู้จะเริ่มตรงไหน หวังว่าคงมีประโยชน์กันบ้าง

Stan จัดว่าเป็น Probabilistic Programming Language อันหนึ่งที่ช่วยให้เราทำโมเดลที่เกี่ยวข้องกับความน่าจะเป็นในแบบที่ต้องอาศัยทฤษฎีเบย์ได้ง่ายขึ้น และปัญหาที่เราสามารถเอา Stan มาใช้ประโยชน์ก็เช่นด้าน optimization ครับ ซึ่งโปรแกรมด้านนี้มันก็มีหลายตัวครับเช่น WinBUGS/OpenBUGS, JAGS และอื่นๆอีกมาก ลองดูเพิ่มเติมที่ http://probabilistic-programming.org/wiki/Home ครับ

ตัวอย่างที่เอามาโชว์นี้ก็มาจากที่ผมเคยพูดไว้แล้วที่  http://www.sakngoi.com/2017/03/22/probabilistic-programming-language/ เรื่องที่เราจะหาค่า \alpha, \beta จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง y = \alpha+\beta x

อันนี้ตัวอย่างcode R สำหรับสร้างข้อมูลที่เอามาfitกับ Stan ครับ

set.seed(12345)
fakedata<-function(alpha, beta, t){
 err<-rnorm(length(t),mean = 0, sd = 5)
 fake<-data.frame(time=t,y=alpha + beta*t +err)
 return(fake)
}
ob<-fakedata(alpha = 1.5, beta = 0.5, 1:100)
plot(ob)

# using lm to fit the fake data
lmfit<-lm(ob$y ~ ob$time)
cat(paste0("alpha = ",format(round(lmfit$coefficients[1],4),nsmall = 4), " beta = ", format(round(lmfit$coefficients[2],4),nsmall = 4)))
abline(lmfit, col='red')

ซึ่งจากข้อมูลที่สร้างมานี้เรามารถใช้คำสั่ง lm หาได้ว่าค่า \alpha, \beta เป็นเท่าไหร่โดยคำตอบที่ได้นี้มาจากเทคนิคที่เรียกว่า QR decomposition ซึ่งจะต่างจากเทคนิค Markov chain Monte Carlo  ที่ใช้ใน Stan ที่คำตอบจะออกมาเป็นในลักษณะของdistributionที่เราประมาณได้ว่าค่าที่เป็นไปได้อาจอยู่แถวค่า median ของdistributionที่ได้นั้น

อันนี้เป็นตัวอย่างcode ที่ผมใช้fitกับข้อมูลด้านบนครับ

library(rstan)
#เป็นfunctionที่สร้างสำหรับให้ค่าเริ่มต้นของตัวแปร
initf<-function(){
 list(alpha=0.1,beta=0.1,sig=0.1)
}
#สร้างตัวแปรdataใหม่ในรูปแบบของlist
dat <- list(y=ob$y,x=ob$time,np=length(ob$y))

#codeของstan
mymodel <- "
functions{
 real fn(real x, real alpha, real beta){
 return alpha + beta*x;
 }
}
data{
 int np;
 vector[np] y;
 vector[np] x;
}
parameters{
 real alpha; 
 real beta;
 real sig;
}
model{
 vector[np] md;

alpha ~ uniform(-5,5);
 beta ~ uniform(-5,5);
 sig ~ uniform(0,30);

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

y ~ normal(md,sig);
}
generated quantities{
 vector[np] y_pred;
 vector[np] md;

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

for(i in 1:np)
 y_pred[i] <- normal_rng(md[i], sig);

}
"
stanout <- stan(model_code = mymodel, data = dat, chains = 1, iter = 10000, init = initf)
summary(stanout,pars=c("alpha","beta"),probs=c(0.025,0.5,0.975))$summary

ตัวอย่างของoutputของค่าที่ได้ครับ

บางคนพอเห็นcodeของ Stanแล้วบอกว่าจะมาเสียเวลาเขียนทำไมยาวๆใช้ lm สิบันทัดเดียวจบเลย ก็จริงครับสำหรับปัญหาง่ายๆแบบนี้ แต่ถ้าปัญหาที่มันซับซ้อนกว่านี้ล่ะ ผมหมายถึงโมเดลที่มันมีหลายตัวแปรที่ซับซ้อนในการคิดและแน่นอนไม่ได้เป็นแบบ linear ล่ะ Stan มันถูกออกแบบมาให้ใช้กับปัญหาที่หลากหลายครับ แต่มันก็ต้องเปลี่ยนวิธีมองปัญหาให้เป็นแบบเบย์ (Bayesian) ครับ อย่างเช่นในปัญหานี้ ….

เดี๋ยวมาต่อ

 

Stan’s Program Blocks

stanblocks

หลักการของ Stan คือจะใช้ตัวแปรหรือ function อะไรก็ต้อง define ก่อน ไม่เหมือนกับ Bugs

บล็อคที่จำเป็นสำหรับ Stan คือ model{} ส่วนอันอื่นเป็นแค่ตัวเสริม แต่ถ้าบล็อคตัวเสริมถูกใส่เข้ามาในโมเดลก็จะต้องเรียงตามลำดับที่แสดงไว้ในรูปคือ เริ่มจาก functions, data, transformed data, parameters, transformed parameters, model และ generated quantities

MCMC ใน R

RStudio-R-JAGS

ถ้าจะทำ Markov Chain Monte Carlo Monte Carlo (MCMC) ใน R นั้นนอกจากจะเขียนแล้วก็ยังมีโปรแกรมช่วยอีกหลายตัวครับ โดยโปรแกรมที่เป็นที่นิยมกันก็ได้แก่ WinBUGS, OpenBUGS, Jags, และก็ Stan ครับ โดย R มี package ที่ช่วยให้เราส่งผ่านหรือรับข้อมูล/โมเดล ระหว่าง R กับโปรแกรม MCMC เหล่านี้ครับ เช่น rjags, RStan, R2WinBUGS, R2OpenBUGS, BRugs ส่วน package ที่ช่วยในการวิเคราะห์ก็อย่างเช่น coda