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
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
Mathematica R Uncategorized

เหตุผลที่ทำไมกินยาartesunateทุก12ชม.ถึงไม่ช่วยเรื่องการรักษามาลาเรีย

เรื่องมันมีอยู่ว่าแบบจำลองเรื่องการดื้อยาartemisininที่เสนอโดย Saralamba et al. (ผมและทีมงาน) ทำนายว่าการเพิ่มขึ้นของระยะเวลาในการกำจัดเชื้อมาลาเรียจากร่างกายคนไข้เป็นผลมาจากการที่ระยะหนึ่งของเชื้อมาลาเรียที่เรียกว่า Rings หรือวงแหวนนั้นมีการตอบสนองกับยาลดลงพร้อมกับเสนอว่าถ้าเปลี่ยนการให้ยาคนไข้ทุก 24 ชม.มาเป็น ทุก12ชม.จะช่วยลดเวลาของการกำจัดเชื้อได้ดีกว่า แต่จากการศึกษาโดย Das et al. ที่ทดลองให้ยาทุก 12 ชม. และเพิ่มขนาดเป็นของยากับคนไข้มาลาเรียพบว่าระยะในการกำจัดเชื้อไม่เปลี่ยนแปลงมากนัก จากผลที่ได้นี้มันแสดงให้เห็นว่าแบบจำลองที่เสนอไปนั้นต้องมีอะไรผิดพลาดแน่นอนและก็เลยทำให้มีการศึกษากันเพิ่มขึ้นพร้อมกับมีการเสนอสมมุติฐานต่างๆเพิ่มขึ้น

ในแบบจำลองอันเดิม เราสมมุติว่าหลังจากที่เชื้อโดนยาแล้วจะตายพร้อมกับถูกกำจัดออกจากร่างการเลย อย่างรวดเร็วโดยม้าม และเป็นเพราะด้วยเหตุผลนี้ที่ทำให้มีการเถียงกันเรื่องการใช้ระยะของการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยาartemisinin ว่าอาจไม่เหมาะสม (https://www.ncbi.nlm.nih.gov/pubmed/26239987) เพราะนั่นหมายความว่าถ้าเชื้อมันยังตอบสนองกับยาได้ดีอยู่ล่ะ แต่ม้ามมันกำจัดเชื้อที่ตายแล้วออกจากร่างการได้ช้ามันก็ดูเหมือนว่าเป็นเชื้อที่ดื้อยาถ้ายังใช้ระยะเวลาในการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยา  แต่มันก็มีหลายงานวิจัยที่แสดงให้เห็นว่าม้ามนั้นกำจัดเชื้อมาลาเรียได้เร็วมากในผู้ป่วยที่รักษาด้วยยาartemisinin (https://www.ncbi.nlm.nih.gov/pubmed/28231817,https://www.ncbi.nlm.nih.gov/pubmed/11992295)

ผมและทีมงานก็เลยเปลี่ยนวิธีที่เชื้อมาลาเรียจะตอบสนองกับยาในแบบจำลองโดยแทนที่เชื้อจะตายเลยหลังจากที่โดนยาก็เปลี่ยนมามาอยู่ในสถานะที่เซลล์เสียหายจากยาแทนแล้วค่อยๆตาย และก็ให้ว่าถ้าเชื้อที่อยู่ในสถานะที่เซลล์เสียหายอยู่ถึงแม้จะโดนยาอีกก็ไม่ได้มีผลเปลี่ยนแปลงอะไรเลย ผลที่ได้จากแบบจำลองก็ดูเหมือนจะสอดคล้องกับทุกผลที่ได้จาก Das et al. นั่นก็คือถึงแม้จะให้ยาทุก 12 ชม.ก็ไม่ได้ลดระยะในการกำจัดเชื้อสักเท่าไหร่ กราฟด้านล่างนี้เป็นตัวอย่างผลที่ได้จากแบบจำลอง สีชมพูแสดงการลดลงของเชื้อในคนไข้ที่กินยาทุก 24 ชม.และสีฟ้าสำหรับคนไข้ที่กินยาทุก12ชม. ซึ่งแสดงให้เห็นว่าไม่ได้แตกต่างกันเท่าไหร่และกราฟนี้ก็แสดงให้เห็นว่าถ้าเราทำการศึกษาแบบที่ Das et al. ทำ โอกาสที่เราจะได้ว่าระยะของการกำจัดเชื้อแบบกินยาทุก24ชม.กับทุก12ชม.จะไม่แตกต่างกันมีสูงถึง 99%

ส่วนด้านล่างนี้ผมทดลองรันแบบจำลองให้ยาวขึ้นเพื่อดูว่ามันจะทำให้ระดับของเชื้อลดลงไปจนถึงเป็นศูนย์หรือไม่โดยสีชมพูคือแบบที่คนไข้กินยาอาทีซูเนตวันล่ะเม็ดและสีฟ้าคือแบบทีกินวันล่ะสองเม็ด จะเห็นได้ว่ามันมีแนวโน้มที่เชื้อจะกลับมาในแบบจำลองนี้ทั้งที่ความเป็นจริงคนไข้เจ้าของข้อมูลนั้นหายจากมาลาเรีย ซึ่งก็เป็นอะไรที่ต้องแก้ไขต่อไป…อาเมน

ในการศึกษานี้ผมใช้ Stan ในการสร้างแบบจำลองและเปรียบเทียบผลที่ได้กับข้อมูลจริงจากคนไข้ครับ

Categories
Uncategorized

ใช้ประโยชน์จาก rstan

ถ้าใครได้ใช้ deSolve สำหรับแก้ปัญหาเชิงตัวเลขของระบบสมการเชิงอนุพันธ์(ode) แล้วจะเห็นว่ามันมีข้อจำกัดเรื่องความเร็วอยู่ เพราะถ้าเขียน function ของสมการด้วย R แล้วจะช้ามาก ยิ่งถ้ามีการเอาfunction ไปใช้ fitting ด้วยเทคนิคอย่าง Markov chain monte carlo แล้วนี่เลิกคิดได้เลยสำหรับระบบใหญ่ ๆ

deSolve เองก็พยายามแก้ปัญหานี้ โดยผู้ใช้สามารถเรียกใช้ function ที่เขียนใน C/C++ ที่คอมไพล์แล้วมาใช้ได้ แต่ก็ยังดูยากพอสมควรสำหรับคนที่ไม่เคยเขียน C/C++ เลย แต่มันก็ยังพอมีวิธีที่ดีกว่าหน่อยที่อยากจะแนะนำสำหรับคนที่อยากจะเพิ่มความเร็วสำหรับใช้กับ deSolve  นั่นก็คือเขียน ในภาษาของ rstan แทนครับ

โดยหลักการก็คือเขียนสมการใน rstan แล้วใช้คำสั่งของ rstan ที่เรียกว่า expose_stan_functions เพื่อที่จะคอมไพล์ function แล้วเราก็จะยังเรียกใช้ได้ใน .GlobalEnv นั่นหมายความว่าเราสามารถเรียกใช้ function ที่คอมไพล์แล้วนี้กับ deSolve ได้ครับ

ตัวอย่างเช่น ผมเขียนสมการ lorenz ในไฟล์ชื่อ lorenz.stan

##lorenz.stan
functions{
  vector lorenz(real t,vector y,vector parms){
    vector[3] ydot;
    // parms[1]=beta ,[2]=sigma,[3]=rho
    
    ydot[1]=parms[2]*(y[2]-y[1]);
    ydot[2]=y[1]*(parms[3]-y[3])-y[2];
    ydot[3]=y[1]*y[2]-parms[1]*y[3];

    return ydot;
  }
}
model{
}

จากนั้นก็ใช้ expose_stan_functions เรียกจากใน R ได้เลยตามนี้ครับ

library(rstan)

expose_stan_functions("lorenz.stan")

# parms[1]=beta ,[2]=sigma,[3]=rho
parameters1<-c(2.66666,10,28)

# x=-10, y=-12, z=30.05
init.state<-c(-10,-12,30.05)

nderivs<-function(t,state,parameters){
  tmp<-lorenz(t,state,parameters)
  return(list(tmp))
}

run1<-ode(y=init.state,times=seq(0,15,0.01), 
   func=nderivs, parms=parameters1)

plot(run1)

 

 

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

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