ถ้าจะทำ Markov Chain Monte Carlo Monte Carlo (MCMC) ใน R นั้นนอกจากจะเขียนแล้วก็ยังมีโปรแกรมช่วยอีกหลายตัวครับ โดยโปรแกรมที่เป็นที่นิยมกันก็ได้แก่ WinBUGS, OpenBUGS, Jags, และก็ Stan ครับ โดย R มี package ที่ช่วยให้เราส่งผ่านหรือรับข้อมูล/โมเดล ระหว่าง R กับโปรแกรม MCMC เหล่านี้ครับ เช่น rjags, RStan, R2WinBUGS, R2OpenBUGS, BRugs ส่วน package ที่ช่วยในการวิเคราะห์ก็อย่างเช่น coda
พอดีเขียนโปรแกรมใช้งาน WinBUGS ผ่าน RLink ใน Mathematica แล้วเจอการส่งและรับค่าแปลกๆของตัวแปรประเภท structure ที่มีการเขียนหรือวางตำแหน่งของอาร์เรย์ไม่เหมือนกัน คือว่า จาก .Data ที่เป็นlistของตัวเลขจะถูกแปลงเป็นอาร์เรย์สองมิติโดยที่ใน Rจะใส่ตัวเลขลงในแต่ละcolumnหรือตามแนวตั้ง แต่WinBUGS จะใส่ตัวเลขตามแถวหรือแนวนอน ระวังกันด้วยเน้อ… http://www.openbugs.net/Manuals/ModelSpecification.html#FormattingOfData
สมมติให้ว่าเรามีผลการทดลองอันหนึ่งที่มีการทำซ้ำอยู่ทั้งหมด 10 ครั้ง และในแต่ละครั้งเราก็จะได้ผลจากการวัดมาทั้งหมด 25 ค่า (t = 1,2,..,25) เมื่อเรานำว่าเขียนกราฟก็จะมีหน้าตาแบบนี้ และเราก็มีว่าโมเดลหรือแบบจำลองที่จะอธิบายผลการทดลองนี้คือ y(t)=C0 exp(-ke*t) ซึ่งถ้าเราทำการ fit ข้อมูลที่ได้จากการทดลองแต่ล่ะครั้งกับโมเดลนี้เราก็จะได้ค่า C0 และ ke สำหรับแต่ละการทดลอง เช่น ถ้าใช้ คำสั่ง NonlinearModelFit หรือ FindFit ก็จะได้ ลองแสดงผลจากการ fit กับข้อมูลจากแต่ละการทดลองด้วย Manipulate หรือจะใช้ WinBUGS ในการหาก็ได้ครับ เช่น ให้ค่าที่ได้จากการวัดจากการทดลอง i ที่เวลา j (y[i,j]) มีการคลาดเคลื่อนจากโมเดล (mod[i,j]) เล็กน้อย (บวก/ลบ SD = 1/sqrt(tau))แบบnormal distribution ลองดู WinBUGS code ด้านล่างนี้นะครับ โดยเรากำหนดให้ว่า ค่า C0 […]
อันนี้เป็นอีกตัวอย่างที่ผมใช้ประโยชน์ 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 ทุกตัว […]
ปัญหาที่เจอคือหลังจากที่ใช้งาน WinBUGS ผ่าน R2WinBUGS ไปสักพักรู้สึกว่าตัว WinBUGS มันช้ากว่าที่ตัวมันเองจะเริ่มประมวนผล ก็พยามยามหาดูว่ามันเป็นเพราะอะไร สุดท้ายไปเจอว่าไฟล์ ชื่อ Registry.odc และ Registry_Rsave.odc ใน System\Rsrc ของ WinBUGS มันมีขนาดใหญ่มากประมาณ 16 GB ! ก็เลยตัดสินใจลบมันทิ้งทั้งคู่ แล้วลอง R2WinBUGS ใหม่ปรากฏว่ามันกลับมาเร็วเหมือนเดิมและเจ้าไฟล์ทั้งสองก็ถูกสร้างขึ้นมาใหม่ที่ขนาดเล็กกว่าเดิม ฉะนั้นใครที่เจอปัญหานี้ให้ลองดูที่ขนาดไฟล์ Registry ทั้งสองนี้ด้วยนะครับ