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
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

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

Categories
Mathematica

ใช้ RLink ช่วยในการทำ MCMC ด้วย WinBUGS

อันนี้เป็นอีกตัวอย่างที่ผมใช้ประโยชน์ RLink ใน Mathematica ครับ

ผมมีแบบจำลองอันหนึ่งที่ผมเขียนด้วยภาษา Oberon (โบราณมาก; http://en.wikipedia.org/wiki/Oberon_(programming_language)) เพื่อใช้ใน WinBUGS อีกทีเพื่อทำ Markov Chain Monte Carlo (MCMC; http://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo) ด้วยเทคนิค Gibbs Sampling (http://en.wikipedia.org/wiki/Gibbs_sampling) ในการประมาณค่าตัวแปรของแบบจำลองครับ

ในการทำ MCMC เพื่อประมาณค่าตัวแปรของแบบจำลองกับข้อมูลแต่ละชุดของผมนี้จะใช้เวลาประมาณ 15 นาที ซึ่งค่อนข้างจะเสียเวลามากถ้าทำทีละครั้งกับข้อมูลทั้งหมดที่มี ดังนั้นผมเลยใช้ RLink ช่วย

วิธีที่ผมทำก็คือใช้ R เรียก WinBUGS ด้วยแพ็คเกจ R2WinBUGS แล้วก็ใช้ RLink ช่วยเรียก R อีกทีครับ จากนั้นผมก็เอามาทำ parallel ใน Mathematica ครับ

อันนี้เป็นตัวอย่าง code ครับ

เราจะทำการโหลด RLink ไปที่ kernels ทุกตัว พร้อมกับเรียกใช้ R2WinBUGS

para1

ด้านล่างนี้เป็น code ที่ผมใช้ในการสร้าง input ไฟล์ของ WinBUGS สำหรับข้อมูลแต่ล่ะชุด

para2

และ อันนี้ code ที่ผมใช้สำหรับ run WinBUGS ครับ

para3

ซึ่งก่อนที่จะ run แบบ parallel นั้นต้องทำให้ทุก ๆ kernel รู้จัก function หรือ ตัวแปรที่จะใช้ก่อนด้วยคำสั่ง DistributeDefinitions

จากนั้นเราก็สามารถที่จะการคำนวณแบบ parallel ได้ครับ

para4

ส่วนอันนี้เป็นตัวอย่างของผลลัพท์ที่ได้จากการทำ MCMC ครับ

para5