แบบจำลองการเพิ่มขึ้นของเชื้อมาลาเรียในคน

กำลังเตรียมdraft งานเก่าที่ทำมาหลายปีแล้วไม่ได้มีโอกาสส่งตีพิมพ์สักที เลยเอามาปัดฝุ่นใหม่โดยใช้ Mathematica รู้สึกได้เลยว่า code มันเขียนไม่กี่บรรทัดเอง แถมแชร์ก็ง่าย https://www.wolframcloud.com/obj/sompob/Published/NJW_Im4a.nb

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

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

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

“All models are wrong, but some are useful”

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

แก้ legend ใน ggplot

สัญลักษณ์จุดกับเส้นของ legend ที่ggplotมันสร้างให้หมุนตาม coordinates ที่flip

วิธีแก้ก็ต้องลงไปแก้ในตอนที่มันสร้าง legend เท่าที่ค้นมาก็จะเขียนประมาณนี้ โดยแก้ที่ตัวแปล $drawkey ของ geom ทั้งหลาย อย่างกรณี geom_point,geom_pointrange ก็จะแก้เป็นเช่น

library(grid)
GeomPointrange$draw_key <-  function (data, params, size)     {
  
  draw_key_vpath <- function (data, params, size) {
    segmentsGrob(0.1, 0.5, 0.9, 0.5, 
                 gp = gpar(lwd = data$size * .pt, lty = data$linetype, 
                           lineend = "butt"), arrow = params$arrow)
  }
  
  grobTree(draw_key_vpath(data, params, size), 
           draw_key_point(transform(data, size = data$size * 4), params))
}

วิธีใช้ก็รัน code นี้ เพื่อแก้ drawkey ก่อนคำสั่ง ggplot

เรียก 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 ตามนี้ครับ