Categories
R Uncategorized

มาลองใช้ Stan แก้ปัญหา

พอดีว่าง ระหว่างรอผลคำนวณอะไรบางอย่าง ผมเห็นมีคนโพสท์ถามที่เวบพันทิปตามนี้ครับ

จากคำถามนี้เราสามารถที่จะใช้ Stan แก้ปัญหาได้ถ้าผมมองว่าการวัดเปรียบเทียบเครื่องมือมาตรฐาน(จริง)กับที่ปรับปรุงขึ้นมานั้นในแต่ล่ะครั้งนั้นไม่ได้เกี่ยวกันเลย และค่าจากวัดของเครื่องที่ปรับปรุงนั้นมีการกระจายรอบค่าจากเครื่องมาตราฐานแบบ normal distribution โดยมีค่า standard deviation หรือ SD อยู่ค่าหนึ่งที่เป็นตัวกำหนดประสิทธิภาพของเครื่อง อย่างเช่น จากข้อมูลที่ให้มาเมื่อวัดเทียบกับเครื่องมือจริงที่วัดได้ 2 แต่ค่าจากเครื่องปรับปรุงวัดมา 3 ครั้งได้ (1.8,2.1,1.9) เราจะสมมุติให้ทั้ง 3 ค่านี้กระจายรอบค่าใดค่าหนึ่ง โดยที่ ถ้าเราปรับปรุงได้เจ๋ง ค่าที่มันกระจายรอบนี้มันก็ควรจะได้เท่ากับค่าจากเครื่องจริงโดยที่มีค่า SD น้อยที่สุดเท่าที่จะเป็นไปได้ แต่ในความเป็นจริงก็จะมีผลผันผวนหรือerrorต่างๆเข้ามาเกี่ยวข้อง ที่ผมจะทำให้ดูนี้เราจะลองใช้โปรแกรมอย่าง Stan มาช่วยหาดูว่าค่ากลางที่เครื่องปรับปรุงหรือสร้างขึ้นมานั้นมันกระจายรอบในแต่ล่ะครั้งของการวัดมันคือค่าอะไรและมี SD เท่าใด

ผมใช้โปรแกรมที่ชื่อ Stan ในภาษา R ผ่าน library ที่เรียกว่า rstan ครับ ด้านล่างนี้ก็เป็น code ที่ผมใช้กับปัญหานี้ครับ

#โหลด library ที่จะใช้งานครับ 
#ในที่นี้มี 2 ตัวคือ rstan กับ bayesplot (ใช้วาดกราฟสรุปผล)
library(rstan)
library(bayesplot)

#เตรียมข้อมูล โดยผมแยกตามค่าของเครื่องจริง
data <- list(
  ob2 = c(1.8,2.1,1.9),
  ob3 = c(3.2,3.0,3.2),
  ob4 = c(3.9,4.3,4.4),
  ob5 = c(5.0, 5.2, 5.5),
  ob6 = c(6.1,5.9,6.5),
  n=3
)

# code ของโมเดล
model <- "
data{
      //จำแหนกประเภทของตัวแปรของข้อมูล ว่าเป็นเลขจำนวนเต็มหรือทศนิยม พร้อมกำหนดขนาด
        int<lower=1> n; 
	real ob2[n];
	real ob3[n];
	real ob4[n];
	real ob5[n];
	real ob6[n];
}
parameters{
// กำหนดประเภทของตัวแปรที่จะใช้ในโมเดล ซึ่งในที่นี้คือค่ากลางที่ข้อมูลมันกระจายรอบ ซึ่งแบ่งตามค่าจริง
// sig เป็นค่า SD 
	real ob2mu;
	real ob3mu;
	real ob4mu;
	real ob5mu;
	real ob6mu;
	real<lower=0> sig;
}
model{
//กำหนดว่าค่ากลางมันและSD อยู่ในช่วงไหน จากการกระจายแบบไหน 
//ซึ่งในที่นี้ผมให้มันมาจาก uniform distribution โดยใส่ช่วงที่คิดว่าค่ามันจะอยู่ในนั้น
	sig ~ uniform(0,3);
	ob2mu ~ uniform(0,7);
	ob3mu ~ uniform(0,7);
	ob4mu ~ uniform(0,7);
	ob5mu ~ uniform(0,7);
	ob6mu ~ uniform(0,7);

// กำหนดให้ค่าที่วัดมาในแต่ล่ะค่าจริงกระจายแบบ normal รอบค่ากลางอันหนึ่ง
	ob2 ~ normal(ob2mu, sig);
	ob3 ~ normal(ob3mu, sig);
	ob4 ~ normal(ob4mu, sig);
	ob5 ~ normal(ob5mu, sig);
	ob6 ~ normal(ob6mu, sig);
 
}
"
fit <- stan(model_code=model, data = data, iter = 10000)

ผลลัพท์ที่ได้ก็จะประมาณนี้ครับ

เห็นได้ว่าค่ากลางหรือ median ที่ 50% นั้นค่อนข้างจะใกล้กับค่าจริงคือ เช่นที่ค่าจริง = 2 เครื่องปรับปรุงทำได้ 1.93(95%CI: (1.63,2.24)) หรือที่ 6 เครื่องปรับปรุงทำได้ที่ 6.15 (CI: (5.86,6.48)) และ SD = 0.25 CI: (0.17-0.44)

posterior <- as.matrix(fit)
plot_title <- ggtitle("Posterior distributions",
                      "with medians and 95% intervals")
mcmc_areas(posterior,
           pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ),
           prob = 0.95) + plot_title
ตัวอย่างการกระจายของค่าที่ได้
posterior2 <- extract(fit, inc_warmup = T, permuted = FALSE)

color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior2,  pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ), 
                n_warmup = 5000,
                facet_args = list(nrow = 6, labeller = label_parsed))
p + facet_text(size = 15)

จากกราฟการกระจายที่ได้จะเห็นได้ว่าผมใส่เป็น uniform distribution ที่ช่วงระหว่าง 0-7 เลยแต่ผลลัพท์ที่ได้เป็น normal distribution ครับ

chains ของแต่ล่ะตัวแปร
Categories
Mathematica Uncategorized

จำกัดช่วงของเลขสุ่มใน Stan

ถ้าจะจำกัดช่วงของเลขที่สุ่มจาก distribution ใดๆสามารถที่จะช่วงที่ต้องการได้โดยใช้ truncation function T ตามนี้ อย่างเช่น ถ้าอยากให้ y สุ่มมาจาก normal ที่ค่าอยู่ในช่วง -0.2 ถึง 1.5 สามารถพิมพ์ได้ตามนี้ครับ

y ~ normal(0,1) T[-0.2,1.5];

หรือจะพิมพ์แบบนี้ก็ได้ครับ

y ~ normal(0,1);
if (y < -0.2 || y > 1.5)
   target += negative_infinity();
else
   target += -log_diff_exp(normal_lcdf(1.5|0,1),normal_lcdf(-0.2|0,1));

แต่ถ้าทำในภาษา Wolfram ก็พิมพ์แบบนี้ได้ครับ

DT = TruncatedDistribution[{-0.2, 1.5}, NormalDistribution[0, 1]]

ls = RandomVariate[DT, 100000];
Histogram[ls]
Categories
Uncategorized

ตัวอย่าง code สำหรับsampling แบบขนานใน cmdstan

Mac OS/Linux

for i in {1..4} 
do 
    ./my_model sample random seed=12345 id=$i data file=my_data output file=samples$i.csv & 
done

Windows

for /l %x in (1, 1, 4) do start /b model sample random seed=12345 id=%x data file=my_data output file=samples%x.csv

เอามาจาก cmdstan-guide-2.20.0.pdf

Categories
Uncategorized

ปัญหา 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 นะครับ

Categories
Uncategorized

ติดตั้ง 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 ครับ

Categories
R

ใช้ Stan ในงาน PKPD

แนะนำครับสำหรับผู้สนใจใช้ Stan ในงาน PKPD modelling

 

 

Categories
R

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) ครับ อย่างเช่นในปัญหานี้ ….

เดี๋ยวมาต่อ

 

Categories
Uncategorized

Stan’s Program Blocks

stanblocks

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

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

Categories
Uncategorized

Stan operators

stan_op

Categories
Uncategorized

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