พอดีว่าง ระหว่างรอผลคำนวณอะไรบางอย่าง ผมเห็นมีคนโพสท์ถามที่เวบพันทิปตามนี้ครับ
จากคำถามนี้เราสามารถที่จะใช้ 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 ครับ
