มันจะกระจายยังไงถ้าไม่ทำอะไรเลย

ทำเล่นๆขำ ว่าถ้าไวรัสที่ว่ามันกระจายไปตามจังหวัดต่างๆที่ว่า จำนวนคนติดเชื้อมันพุ่งเยอะขนาดไหน แบบว่าถ้าคนมันคิดว่าเป็นหวัดธรรมดา

โห … ตัวเลขพุ่งขึ้นเร็วมากกก ถ้าเห็นโมเดลทำนายเรื่องอะไรก็ตามอยากให้คิดกันเสมอว่า

“All models are wrong, but some are useful”

ฉะนั้นอย่าเพิ่งไปเชื่ออะไรง่ายๆครับ ทุกโมเดลมันมีสมมุติฐานในการสร้างของมัน ระบบในธรรมชาติมันซับซ้อนมากครับ ไม่มีใครที่จะโมเดลมันได้หมดหรอกครับ ทุกโมเดลมันจะตัดความยุ่งยาก ซับซ้อนออกเพื่อให้โมเดลได้ง่ายครับ

มาลองใช้ 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 ของแต่ล่ะตัวแปร

เรียก compiled function จาก compiled function

โดยปกติ Mathematica หรือภาษา Wolfram จะมีปัญหากับการเรียก function ที่สร้างจาก Compile ใน function ที่จะสร้างจาก Compile อีกทีได้ (ส่วนใหญ่เป็นเรื่องที่เกี่ยวกับว่า Funcions ไหนที่มัน compile ได้หรือไม่ได้) แต่ใน version 12 นี้เราสามารถที่จะทำแบบนั้นได้ ถ้า function เหล่านั้นถูก compile ด้วยคำสั่ง FunctionCompile ซึ่งผมมองว่ามันสะดวกอย่างมาก และที่สำคัญเราสามารถที่จะ export สร้างเป็น library ไฟล์ได้ด้วย โดยคำสั่ง FunctionCompileExport สกุลไฟล์ที่สร้างได้ก็ได้แก่

ตัวอย่างการเรียกใช้ compiled function ใน compiled function อีกที ก็สมมุติให้ผมมี function สำหรับทำ bisection เพื่อหาค่าที่ทำให้ function ที่สนใจเป็น 0 หรือใกล้ 0 มากๆ

BisectionMethod = FunctionCompile[
    Function[
        {Typed[f, {"Real64"} -> "Real64"], Typed[lim0, "Real64"], Typed[lim1, "Real64"]},
      Module[{Maxiter = 100, tol = 10.^-8, iter = 0, x0 = lim0, x1 = lim1, mid = 0.0, f0, f1, fmid},
      f0 = f[x0]; 
      f1 = f[x1];
            While[Abs[x1 - x0] >= tol && iter++ < Maxiter,
                mid = (x1 + x0)/2.;
                fmid = f[mid];
                If[Xor[fmid > 0, f0 > 0],
                    x1 = mid; f1 = fmid,
                    x0 = mid; f0 = fmid]];
               mid]]]

โดยfunction ที่ผมจะหา root คือ f(x) = sin(x)+exp(x)

f = FunctionCompile[Function[Typed[arg, "Real64"], Sin[arg] + Exp[arg]]]

จะเห็นได้ว่าทั้ง BisectionMethod และ f ถูกสร้างจาก FunctionCompile แต่ f จะถูกเรียกเข้าไปใช้ใน BisectionMethod อีกที

BisectionMethod[f, -1.0, 1.0]-FindRoot[Sin[x] + Exp[x], {x, -1}][[1, 2]]
(* -5.18667*10^-11 *)

การแพร่กระจายของโรคไข้เลือดออกในจีน

พอดีว่าต้องทำ slides ไปโชว์ว่างานที่รับอยู่ไปถึงไหนแล้ว มันมีส่วนหนึ่งที่ต้องทำให้เห็นว่าผลที่ได้จากโมเดลการระบาดของไข้เลือดออกในจีนกับข้อมูลจริงนั้นมันใกล้เคียงกัน พอดีผมไปเจองานอยู่งานหนึ่งที่น่าสนใจคือ https://bmcmedicine.biomedcentral.com/articles/10.1186/s12916-015-0336-1

Figure 5

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

โดยผมเขียน code ตามนี้ครับ

เรื่องการเรียกใช้ package ใน package อีกทีของภาษา Wolfram

พอดีกำลังรับงานพัฒนาโมเดลตัวหนึ่งด้วยภาษา Wolfram แล้วเจอปัญหาเรื่องการเรียกใช้ Package ในอีก Package หนึ่ง code ที่เขาส่งมาให้ช่วยดูให้มันมีการเรียกใช้ package ตั้งแต่คำสั่ง BeginPackage เช่น

BeginPackage[ "Package`", {"Package1`", "Package2`"}]
  ...

แต่ปรากฏว่าหลังจากที่ทำการ deploy ไปแล้ว เจอว่ามันมีการคำนวณอะไรบางอย่างผิด ทั้งที่ในตัว unit test ก็ไม่เจอ ผมก็นั่งไล่ดูทีล่ะบรรทัดจนพบว่ามันมีการเขียนทับ functions กัน เพราะชื่อดันเหมือนกัน

ผมก็เลยเปลี่ยนวิธีการเรียก packages ด้วย Needs ไปใช้วิธีเรียกทีล่ะตัวหลังคำสั่ง BeginPackage แทน เช่น

BeginPackage[ "Package`"]
  Needs[ "Package1`"]
  Needs[ "Package2`"]
  ...

ด้วยวิธีนี้ตัว คำสั่งจาก package1 หรือ package2 จะถูกเรียกใช้ได้ภายใน package แต่ด้านนอก package จะเรียกไม่ได้ พูดง่ายๆคือไม่มีผลอะไรกับตัวแปรนอก package