จากที่เคยไปโม้ไว้เยอะว่า 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/ เรื่องที่เราจะหาค่า 
จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง 
อันนี้ตัวอย่าง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 หาได้ว่าค่า เป็นเท่าไหร่โดยคำตอบที่ได้นี้มาจากเทคนิคที่เรียกว่า 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) ครับ อย่างเช่นในปัญหานี้ ….
เดี๋ยวมาต่อ
Like this:
Like Loading...
จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง 
