SakNgoi

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

จากคำถามนี้เราสามารถที่จะใช้ 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 ครับ