เรื่องปวดหัว R structure vs WinBUGS structure

พอดีเขียนโปรแกรมใช้งาน WinBUGS  ผ่าน RLink ใน Mathematica แล้วเจอการส่งและรับค่าแปลกๆของตัวแปรประเภท structure ที่มีการเขียนหรือวางตำแหน่งของอาร์เรย์ไม่เหมือนกัน คือว่า จาก .Data ที่เป็นlistของตัวเลขจะถูกแปลงเป็นอาร์เรย์สองมิติโดยที่ใน Rจะใส่ตัวเลขลงในแต่ละcolumnหรือตามแนวตั้ง แต่WinBUGS จะใส่ตัวเลขตามแถวหรือแนวนอน

R-structurewb-structureระวังกันด้วยเน้อ…

http://www.openbugs.net/Manuals/ModelSpecification.html#FormattingOfData

 

ใช้ Mathematica อ่านไฟล์ coda ที่ได้จาก WinBUGS/OpenBUGS

เราใช้Mathematicaอ่านไฟล์ coda กับ index ของมันที่ได้จาก WinBUGS/OpenBUGS แต่ล่ะ chain ด้วยคำสั่ง ReadList เช่น
codaIndex=Partition[ReadList[“codaIndex.txt”,Word],3];
codaCH1=Partition[ReadList[“coda1.txt”,Word],2];

จากนั้นก็เขียนคำสั่งตามนี้เลยครับ โดยที่คำสั่ง posCodaIdx จะเป็นการหาตำแหน่งของโหนด(obj)ที่เราต้องการอยู่ตำแหน่งที่เท่าใดถึงเท่าใดบ้างในไฟล์coda ส่วน codaObj จะเป็นการแสดงลิสต์ตัวเลขหรือchainของobjนั้น

posCodaIdx[obj_]:=ToExpression@(codaIndex[[Position[codaIndex,obj]//Flatten//First]]//Rest)

codaObj[obj_]:=ToExpression@codaCH1[[#1;;#2,2]]&@@posCodaIdx[obj]

ตัวอย่าง

codaex

 

การ fit ข้อมูลหลายชุดกับโมเดล

สมมติให้ว่าเรามีผลการทดลองอันหนึ่งที่มีการทำซ้ำอยู่ทั้งหมด 10 ครั้ง และในแต่ละครั้งเราก็จะได้ผลจากการวัดมาทั้งหมด 25 ค่า (t = 1,2,..,25) เมื่อเรานำว่าเขียนกราฟก็จะมีหน้าตาแบบนี้

noo1

และเราก็มีว่าโมเดลหรือแบบจำลองที่จะอธิบายผลการทดลองนี้คือ
y(t)=C0 exp(-ke*t)
ซึ่งถ้าเราทำการ fit ข้อมูลที่ได้จากการทดลองแต่ล่ะครั้งกับโมเดลนี้เราก็จะได้ค่า C0 และ ke สำหรับแต่ละการทดลอง เช่น ถ้าใช้ คำสั่ง NonlinearModelFit หรือ FindFit ก็จะได้

noo2

ลองแสดงผลจากการ fit กับข้อมูลจากแต่ละการทดลองด้วย Manipulate

noo3

หรือจะใช้ WinBUGS ในการหาก็ได้ครับ เช่น ให้ค่าที่ได้จากการวัดจากการทดลอง i ที่เวลา j (y[i,j]) มีการคลาดเคลื่อนจากโมเดล (mod[i,j]) เล็กน้อย (บวก/ลบ SD = 1/sqrt(tau))แบบnormal distribution
ลองดู WinBUGS code ด้านล่างนี้นะครับ โดยเรากำหนดให้ว่า ค่า C0 กับ ke ถูกสุ่มมาจาก uniform distribution ตาม range ที่กำหนดตาม code

noo4

ถ้าใครจะใช้ Mathematica เรียก WinBUGS อีกทีก็สามารถทำได้ตามนี้ครับ

noo5

โดยผลที่ได้ก็จะประมาณนี้ครับ

noo6

ใครที่สนใจเกี่ยวกับวิธีการที่ WinBUGS ใช้นั้นลองดูเรื่อง Bayesian Statistics , Markov Chain Monte Carlo (MCMC) ครับ

ซึ่งในการหาค่า C0 กับ ke แบบนี้เรามองว่าแต่ละการทดลองไม่มีขึ้นต่อกันหรือไม่สนใจว่า C0 หรือ ke จากแต่ละการทดลองมีความสัมพันธ์กันอย่างไร แต่ถ้าเรามองว่า C0 หรือ ke จากการทดลองทั้งหมดนี้มีความสัมพันธ์กันอย่างเช่น C0 หรือ ke จากแต่ละการทดลองอาจจะเป็นค่าที่สุ่มมาจาก normal distribution อันหนึ่งที่มีค่า mean หรือ sd อันหนึ่ง เราก็สามารถหาได้เช่น

noo7

โดยผลที่ได้ก็ตามนี้ครับ

noo8

ซึ่งวิธีคิดแบบข้างบนนี้เขาจะเรียกว่า multilevel หรือ hierachical คือเราสนใจว่าในระดับ population ของข้อมูล เจ้าตัวพารามิเตอร์ที่สนใจนั้นมันฟอร์มกันเป็น distribution อะไรหรือเปล่า หรือสุ่มมาจาก distribution อันใดอันหนึ่งหรือไม่ ซึ่งเราก็ต้องออกแบบโมเดลของเราตามข้อมูล หรือความสัมพันธ์ที่เรารู้ ถ้าสนใจเรื่องนี้ผมแนะนำหนังสือเล่มนี้ครับ

Data Analysis Using Regression and Multilevel/Hierarchical Models โดย Andrew Gelman และ Jennifer Hill

เดี๋ยวมาต่อ..

ใช้ 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

 

ปัญหา R2WinBUGS เรียก WinBUGS ช้า

ปัญหาที่เจอคือหลังจากที่ใช้งาน WinBUGS ผ่าน R2WinBUGS ไปสักพักรู้สึกว่าตัว WinBUGS มันช้ากว่าที่ตัวมันเองจะเริ่มประมวนผล ก็พยามยามหาดูว่ามันเป็นเพราะอะไร สุดท้ายไปเจอว่าไฟล์ ชื่อ Registry.odc และ Registry_Rsave.odc ใน System\Rsrc ของ WinBUGS มันมีขนาดใหญ่มากประมาณ 16 GB ! ก็เลยตัดสินใจลบมันทิ้งทั้งคู่ แล้วลอง R2WinBUGS ใหม่ปรากฏว่ามันกลับมาเร็วเหมือนเดิมและเจ้าไฟล์ทั้งสองก็ถูกสร้างขึ้นมาใหม่ที่ขนาดเล็กกว่าเดิม ฉะนั้นใครที่เจอปัญหานี้ให้ลองดูที่ขนาดไฟล์ Registry ทั้งสองนี้ด้วยนะครับ

การเพิ่ม function เข้าไปใน WINBUGS

การเพิ่ม function เข้าไปใน WINBUGS

ผมมีข้อมูลที่ต้องการนำมาเปรียบเทียบกับแบบจำลองการเพิ่มจำนวนของเชื้อมาลาเรียในร่างกายผู้ติดเชื้อที่เสนอโดย N.J. White et al. ด้วยเทคนิค Markov chain Monte Carlo (MCMC) และผมได้เขียนโปรแกรมด้วย Mathematica ตัวอย่างตามนี้ครับ

เพื่อช่วยในการเปรียบเทียบนี้แต่ติดปัญหาว่าข้อมูลมีจำนวนมากและโปรแกรมที่ผมเขียนด้วย Mathematica นี้ค่อนข้างช้า มีคนแนะนำให้ทดลองใช้ WINBUGS แต่ภาษาของ BUGS เองไม่เหมาะที่จะนำมาประยุกต์ใช้กับแบบจำลองที่ผมต้องการได้  แต่ผมก็ไปพบว่าเราสามารถที่จะเพิ่มคำสั่งหรือfunction เข้าไปในตัว WINBUGS ได้ด้วยการเขียนภาษาที่เรียกว่า Component Pascal  ผ่านโมดูลอันหนึ่งที่เรียกว่า WBDev ด้วยตัวโปรแกรมที่เรียกว่า BlackBox  ซึ่งผมจะอธิบายคร่าวๆว่าทำอย่างไร
1. ก่อนอื่นเลยก็ทำการติดตั้ง BlackBox จากนั้นก็ copy ไฟล์และ folder ทั้งหมดของ WINBUGS (version > 1.4) ไปว่างไว้ใน folder ของ BlackBox
2. จากนั้นก็เปิด template ไฟล์ที่ชื่อ VectorTemplate.odc หรือ ScalarTemplate.odc ใน BlackBox จาก BlackBoxFolder\WBDev\Mod ขึ้นกับว่า function ที่เขียนมีการ return ค่ากลับเป็นเลขโดด ๆ (scalar) หรือเป็นอาร์เรย์ของตัวเลข (vector)
3. ทำการแก้ไข template ตามต้องการ จากนั้นก็บันทึกไฟล์ไว้ที่เดียวกันกับ template ที่  BlackBoxFolder\WBDev\Mod (อย่าลืมเปลี่ยนชื่อนะครับ)
4. ให้เปิดไฟล์ชื่อ Functions.odc จาก BlackBoxFolder\WBDev\Rsrc แล้วเพิ่มชื่อ function ของเราเข้าไป ดูตัวอย่างในไฟล์นั้นได้ครับ
5. ปิด BlackBox แล้วเปิดใหม่ครับ จะเห็นว่าเราสามารถใช้ WINBUGS จากใน BlackBox ได้เลยพร้อมกัน function ที่เราเพิ่มเข้าไปครับ
ขอให้สนุกกับ MCMC ครับ