SakNgoi

Tag: RStan

จำกัดช่วงของเลขสุ่มใน Stan

ถ้าจะจำกัดช่วงของเลขที่สุ่มจาก distribution ใดๆสามารถที่จะช่วงที่ต้องการได้โดยใช้ truncation function T ตามนี้ อย่างเช่น ถ้าอยากให้ y สุ่มมาจาก normal ที่ค่าอยู่ในช่วง -0.2 ถึง 1.5 สามารถพิมพ์ได้ตามนี้ครับ หรือจะพิมพ์แบบนี้ก็ได้ครับ แต่ถ้าทำในภาษา Wolfram ก็พิมพ์แบบนี้ได้ครับ

ปัญหา Path ของ Rtools

rtools มันจำเป็นสำหรับ r package ที่มี code จากภาษาอื่น อย่าง c/c++ หรือ fortran ในที่จะต้องมีการ compile ระหว่างการติดตั้ง โดยปกติแล้วมัน จะถูกลงไว้ที่ PATH=”C:\Rtools\bin;${PATH}” ถ้าลงไว้ที่อื่นและต้องการเรียกใช้ใน R ตัวแปรชื่อ BINPREF จะเป็นตัวแปรที่ R จะเรียกหา path ของ rtools ฉะนั้นถ้าจะเรียกใช้ rtools ที่ลงไว้ที่อื่นก็ต้องเซ็ต path ให้กับตัวแปรนี้ เช่น BINPREF=”X:/R/Rtools-3.5/mingw_$(WIN)/bin/” ** สังเกตว่า BINPREF ใช้ forward slash นะครับ

Stan กับ ปัญหา linear regression

จากที่เคยไปโม้ไว้เยอะว่า Stan มันเจ๋งยังไงให้กับเพื่อนร่วมงานฟัง เมื่อวันศุกร์ที่แล้วผมก็เลยนัดโชว์แบบคร่าวๆพร้อมกับสาธิตการใช้งานกับปัญหาง่ายๆอย่าง linear regression ให้เค้าดูกันครับ ก็เลยคิดว่าจะเอาที่ไปโม้ไว้มาใส่ไว้ที่นี่ด้วยเผื่อว่าใครสนใจใช้งาน Stan แต่ไม่รู้จะเริ่มตรงไหน หวังว่าคงมีประโยชน์กันบ้าง Stan จัดว่าเป็น Probabilistic Programming Language อันหนึ่งที่ช่วยให้เราทำโมเดลที่เกี่ยวข้องกับความน่าจะเป็นในแบบที่ต้องอาศัยทฤษฎีเบย์ได้ง่ายขึ้น และปัญหาที่เราสามารถเอา Stan มาใช้ประโยชน์ก็เช่นด้าน optimization ครับ ซึ่งโปรแกรมด้านนี้มันก็มีหลายตัวครับเช่น WinBUGS/OpenBUGS, JAGS และอื่นๆอีกมาก ลองดูเพิ่มเติมที่ http://probabilistic-programming.org/wiki/Home ครับ ตัวอย่างที่เอามาโชว์นี้ก็มาจากที่ผมเคยพูดไว้แล้วที่  http://www.sakngoi.com/2017/03/22/probabilistic-programming-language/ เรื่องที่เราจะหาค่า จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง อันนี้ตัวอย่างcode R สำหรับสร้างข้อมูลที่เอามาfitกับ Stan ครับ set.seed(12345) fakedata<-function(alpha, beta, t){ err<-rnorm(length(t),mean = 0, sd = 5) fake<-data.frame(time=t,y=alpha + beta*t +err) return(fake) } ob<-fakedata(alpha = 1.5, beta = 0.5, […]

เหตุผลที่ทำไมกินยาartesunateทุก12ชม.ถึงไม่ช่วยเรื่องการรักษามาลาเรีย

เรื่องมันมีอยู่ว่าแบบจำลองเรื่องการดื้อยาartemisininที่เสนอโดย Saralamba et al. (ผมและทีมงาน) ทำนายว่าการเพิ่มขึ้นของระยะเวลาในการกำจัดเชื้อมาลาเรียจากร่างกายคนไข้เป็นผลมาจากการที่ระยะหนึ่งของเชื้อมาลาเรียที่เรียกว่า Rings หรือวงแหวนนั้นมีการตอบสนองกับยาลดลงพร้อมกับเสนอว่าถ้าเปลี่ยนการให้ยาคนไข้ทุก 24 ชม.มาเป็น ทุก12ชม.จะช่วยลดเวลาของการกำจัดเชื้อได้ดีกว่า แต่จากการศึกษาโดย Das et al. ที่ทดลองให้ยาทุก 12 ชม. และเพิ่มขนาดเป็นของยากับคนไข้มาลาเรียพบว่าระยะในการกำจัดเชื้อไม่เปลี่ยนแปลงมากนัก จากผลที่ได้นี้มันแสดงให้เห็นว่าแบบจำลองที่เสนอไปนั้นต้องมีอะไรผิดพลาดแน่นอนและก็เลยทำให้มีการศึกษากันเพิ่มขึ้นพร้อมกับมีการเสนอสมมุติฐานต่างๆเพิ่มขึ้น ในแบบจำลองอันเดิม เราสมมุติว่าหลังจากที่เชื้อโดนยาแล้วจะตายพร้อมกับถูกกำจัดออกจากร่างการเลย อย่างรวดเร็วโดยม้าม และเป็นเพราะด้วยเหตุผลนี้ที่ทำให้มีการเถียงกันเรื่องการใช้ระยะของการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยาartemisinin ว่าอาจไม่เหมาะสม (https://www.ncbi.nlm.nih.gov/pubmed/26239987) เพราะนั่นหมายความว่าถ้าเชื้อมันยังตอบสนองกับยาได้ดีอยู่ล่ะ แต่ม้ามมันกำจัดเชื้อที่ตายแล้วออกจากร่างการได้ช้ามันก็ดูเหมือนว่าเป็นเชื้อที่ดื้อยาถ้ายังใช้ระยะเวลาในการกำจัดเชื้อเป็นตัวบอกถึงการดื้อยา  แต่มันก็มีหลายงานวิจัยที่แสดงให้เห็นว่าม้ามนั้นกำจัดเชื้อมาลาเรียได้เร็วมากในผู้ป่วยที่รักษาด้วยยาartemisinin (https://www.ncbi.nlm.nih.gov/pubmed/28231817,https://www.ncbi.nlm.nih.gov/pubmed/11992295) ผมและทีมงานก็เลยเปลี่ยนวิธีที่เชื้อมาลาเรียจะตอบสนองกับยาในแบบจำลองโดยแทนที่เชื้อจะตายเลยหลังจากที่โดนยาก็เปลี่ยนมามาอยู่ในสถานะที่เซลล์เสียหายจากยาแทนแล้วค่อยๆตาย และก็ให้ว่าถ้าเชื้อที่อยู่ในสถานะที่เซลล์เสียหายอยู่ถึงแม้จะโดนยาอีกก็ไม่ได้มีผลเปลี่ยนแปลงอะไรเลย ผลที่ได้จากแบบจำลองก็ดูเหมือนจะสอดคล้องกับทุกผลที่ได้จาก Das et al. นั่นก็คือถึงแม้จะให้ยาทุก 12 ชม.ก็ไม่ได้ลดระยะในการกำจัดเชื้อสักเท่าไหร่ กราฟด้านล่างนี้เป็นตัวอย่างผลที่ได้จากแบบจำลอง สีชมพูแสดงการลดลงของเชื้อในคนไข้ที่กินยาทุก 24 ชม.และสีฟ้าสำหรับคนไข้ที่กินยาทุก12ชม. ซึ่งแสดงให้เห็นว่าไม่ได้แตกต่างกันเท่าไหร่และกราฟนี้ก็แสดงให้เห็นว่าถ้าเราทำการศึกษาแบบที่ Das et al. ทำ โอกาสที่เราจะได้ว่าระยะของการกำจัดเชื้อแบบกินยาทุก24ชม.กับทุก12ชม.จะไม่แตกต่างกันมีสูงถึง 99% ส่วนด้านล่างนี้ผมทดลองรันแบบจำลองให้ยาวขึ้นเพื่อดูว่ามันจะทำให้ระดับของเชื้อลดลงไปจนถึงเป็นศูนย์หรือไม่โดยสีชมพูคือแบบที่คนไข้กินยาอาทีซูเนตวันล่ะเม็ดและสีฟ้าคือแบบทีกินวันล่ะสองเม็ด จะเห็นได้ว่ามันมีแนวโน้มที่เชื้อจะกลับมาในแบบจำลองนี้ทั้งที่ความเป็นจริงคนไข้เจ้าของข้อมูลนั้นหายจากมาลาเรีย ซึ่งก็เป็นอะไรที่ต้องแก้ไขต่อไป…อาเมน ในการศึกษานี้ผมใช้ Stan […]

ใช้ประโยชน์จาก rstan

ถ้าใครได้ใช้ deSolve สำหรับแก้ปัญหาเชิงตัวเลขของระบบสมการเชิงอนุพันธ์(ode) แล้วจะเห็นว่ามันมีข้อจำกัดเรื่องความเร็วอยู่ เพราะถ้าเขียน function ของสมการด้วย R แล้วจะช้ามาก ยิ่งถ้ามีการเอาfunction ไปใช้ fitting ด้วยเทคนิคอย่าง Markov chain monte carlo แล้วนี่เลิกคิดได้เลยสำหรับระบบใหญ่ ๆ deSolve เองก็พยายามแก้ปัญหานี้ โดยผู้ใช้สามารถเรียกใช้ function ที่เขียนใน C/C++ ที่คอมไพล์แล้วมาใช้ได้ แต่ก็ยังดูยากพอสมควรสำหรับคนที่ไม่เคยเขียน C/C++ เลย แต่มันก็ยังพอมีวิธีที่ดีกว่าหน่อยที่อยากจะแนะนำสำหรับคนที่อยากจะเพิ่มความเร็วสำหรับใช้กับ deSolve  นั่นก็คือเขียน ในภาษาของ rstan แทนครับ โดยหลักการก็คือเขียนสมการใน rstan แล้วใช้คำสั่งของ rstan ที่เรียกว่า expose_stan_functions เพื่อที่จะคอมไพล์ function แล้วเราก็จะยังเรียกใช้ได้ใน .GlobalEnv นั่นหมายความว่าเราสามารถเรียกใช้ function ที่คอมไพล์แล้วนี้กับ deSolve ได้ครับ ตัวอย่างเช่น ผมเขียนสมการ lorenz ในไฟล์ชื่อ lorenz.stan […]

Stan’s Program Blocks

หลักการของ Stan คือจะใช้ตัวแปรหรือ function อะไรก็ต้อง define ก่อน ไม่เหมือนกับ Bugs บล็อคที่จำเป็นสำหรับ Stan คือ model{} ส่วนอันอื่นเป็นแค่ตัวเสริม แต่ถ้าบล็อคตัวเสริมถูกใส่เข้ามาในโมเดลก็จะต้องเรียงตามลำดับที่แสดงไว้ในรูปคือ เริ่มจาก functions, data, transformed data, parameters, transformed parameters, model และ generated quantities

MCMC ใน R

ถ้าจะทำ Markov Chain Monte Carlo Monte Carlo (MCMC) ใน R นั้นนอกจากจะเขียนแล้วก็ยังมีโปรแกรมช่วยอีกหลายตัวครับ โดยโปรแกรมที่เป็นที่นิยมกันก็ได้แก่ WinBUGS, OpenBUGS, Jags, และก็ Stan ครับ โดย R มี package ที่ช่วยให้เราส่งผ่านหรือรับข้อมูล/โมเดล ระหว่าง R กับโปรแกรม MCMC เหล่านี้ครับ เช่น rjags, RStan, R2WinBUGS, R2OpenBUGS, BRugs ส่วน package ที่ช่วยในการวิเคราะห์ก็อย่างเช่น coda

Back to top
%d bloggers like this: