Tag: R

จากที่เคยเขียนไว้เกี่ยวกับ r package อันหนึ่งที่ผมเขียนชื่อ maemod สำหรับช่วยให้คนที่สนใจอยากคำนวณพวก ode ได้ง่ายขึ้น (ดูข้างล่าง) maemod (แม่มด) มีคนสนใจว่าถ้าพวกตัวแปร state ต่างๆนั้นเป็นแบบ array จะทำอย่างไร ผมเลยเขียนตัวอย่างพร้อมกับเพิ่มความสามารถด้านarrayนี้เข้าไป พอใช้ได้ไปคร่าวๆก่อน ดูตัวอย่างข้างล่างนี้ครับ # Example from Berkeley Madonna # for using Array # METHOD RK4 # # STARTTIME = 0 # STOPTIME = 20 # DT = 0.02 # # d/dt (A[1]) = -k[1]*A[1] # d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]

จากที่เคยไปโม้ไว้เยอะว่า 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,

ถ้าใครได้ใช้ 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

ตัวอย่างการทำเมฆกลุ่มคำใน R ครับ (Windows10 + R > 3.2 + RStudio 0.99.903) library ที่จะใช้มีสองตัวคือ RLongLexTo (ดูวิธีติดตั้งที่ https://github.com/slphyx/RLongLexTo) กับ wordcloud2 (https://github.com/Lchiffon/wordcloud2) library(RLongLexTo) library(wordcloud2) เซ็ต LC_CTYPE ให้ใช้ภาษาไทยครับ Sys.setlocale("LC_CTYPE","Thai") จากนั้นก็สร้างtext file ของข้อความที่ต้องการสร้างเมฆกลุ่มคำ ในที่นี้ผมcopyข้อความมาจาก http://king.kapook.com/royal_words_2545.php แล้วสร้างเป็นไฟล์ชื่อtest.txt อันนี้ทำง่ายๆเลยครับโดยการลากเม้าส์ ไฮไลท์ข้อความที่ต้องการบนweb browser (MS Edge) แล้วกด ctlr+c เพื่อcopy แล้วก็เอาไปวางในnotepad (ctrl+v) แล้วก็save โดยที่ตอนเซฟผมเลือก Encoding เป็น ANSI จากนั้นก็ใช้คำสั่ง RLongLexToF เพื่อทำการแยกคำโดยคำสั่งนี้จะสร้างไฟล์outputออกมา(ผมให้ชื่อ outtest.txt ครับ) โดยที่แต่ล่ะคำจะแยกกันโดยมีเครื่องหมาย “|” คั่นอยู่ RLongLexToF(inputfilename = "test.txt",outputfilename = "outtest.txt") จากนั้นก็อ่านไฟล์ที่แยกคำแล้วมาสร้างเป็นvector ตามนี้ครับ outtxt<-as.vector(strsplit(readLines("outtest.txt"),"[|]")[[1]])

กำลังเขียน library ที่แปลงสมการ ode เป็นfunctionใน R ปรากฏว่าระหว่างที่ทดสอบสมการซึ่งยาวพอสมควร(ประมาณ 40 สมการ) function ที่เขียนไม่สามารถทำตามที่สั่งได้ มีบางส่วนของสมการถูกตัดไป เลยค้นข้อมูลเจอว่า console ของ R จะ limit ความยาวของคำสั่งไว้ที่ 4095 bytes (ไม่ใช่จำนวนตัวอักษร) https://cran.r-project.org/doc/manuals/R-intro.html#R-commands_003b-case-sensitivity-etc

ก่อนติดตั้ง package PIO https://github.com/rietho/IPO/blob/master/vignettes/IPO.Rmd นี้ ใน R ให้ yum install netcdf netcdf-devel ก่อน ครับแล้วชีวิตจะง่ายขึ้น 🙂

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

จาก http://wp.me/p1rsA6-pn ผมเอา code ที่เขียนมาทำเป็น r-package แล้วเรียกมันว่า RLongLexTo (https://github.com/slphyx/RLongLexTo) ครับเผื่อว่าใครอยากเอาไปพัฒนาต่อ  เบื้องต้นนี้ผมทดสอบเฉพาะสำหรับใช้งานบนวินโดว์ 10 เท่านั้น วิธีติดตั้ง *RLongLexTo ต้องการบางคำสั่งจาก package rJava ครับ ดังนั้นผู้ใช้ก็ต้องลง rJava ก่อนครับ library(devtools) install_github("slphyx/RLongLexTo") ##คำแนะนำสำหรับผู้ที่ลง rJava แล้วมีปัญหาใช้RLongLexToไม่ได้ ผมแนะนำให้ลง java JDK (http://www.oracle.com/technetwork/java/javase/downloads/index.html) หรือไม่ก็ลองดูเรื่อง PATH ของ Java ก่อนครับ ใน RLongLexTo นี้จะมีคำสั่งหลักคือ RLongLexToC ที่จะทำการรับข้อความภาษาไทยมาเพื่อทำการแบ่งคำ วิธีการใช้งาน RLongLexToC("ข้อความที่ต้องการแบ่งคำ") หากมีปัญหาเรื่องการแสดงผลภาษาไทยก็ให้ลองเซ็ตlocaleเป็นภาษาไทยก่อนใช้งานอีกทีครับ Sys.setlocale(locale="Thai") มีหลายคนเมล์มาสอบถามเรื่องการเอา RLongLexTo ไปใช้พร้อมกับถามเรื่องรายละเอียดของอัลกอรึทึ่มด้วย ซึ่งผมก็อยากบอกตรงนี้เลยนะครับว่าผมไม่รู้เรื่องหรอกครับ  ผมเพียงแต่เอา code ต้นฉบับมาดัดแปลงให้ใช้ใน R ตามที่มีคนถามมาและผมก็ว่างทำให้ดูเป็นตัวอย่างเท่านั้นครับ ส่วนตัวแล้วผมอยากจะบอกว่าคำสั่งสำหรับจัดการเรื่องข้อมูลประเภท string จาก package ชื่อ stringr

มีคนส่งข้อความมาว่าให้ช่วยทำให้ดูหน่อยว่าจะเอาโปรแกรมตัดคำภาษาไทย http://www.sansarn.com/lexto/ ซึ่งเขียนด้วยภาษา Java โดย Choochart Haruechaiyasak มาใช้ใน R ได้อย่างไร ด้วยความที่เห็นว่าผู้ถามเป็นหญิง แถมส่งการ์ตูนมาด้วยว่าร้องไห้หนักมาก ไอ้เราก็เป็นคนใจอ่อนด้วยสิ เลยทำให้ดูแบบคร่าวๆตามนี้ครับ R สามารถที่เรียกโปรแกรมที่เขียนด้วย Java มาใช้ใน R หมายถึงมีการส่งผ่านค่าหรือรับค่าต่างๆจากโปรแกรมนั้นได้โดยใช้คำสั่งที่อยู่ในpackageอย่าง rJava ครับรายละเอียดก็ไปหาอ่านเอาเองนะครับ โปรแกรม LexTo ที่ให้โหลดกันนั้นจะมีfunction main อยู่ในไฟล์ LongLexTo.java ครับ ผมก็แก้ไขอะไรบางอย่าง(จำไม่ได้ ต้องเปิดดู) หลักๆก็วิธีที่เขียนผลลัพท์ออกมา แล้วก็คอมไพล์ใหม่ซึ่งเราก็จะได้เป็น .class มาครับ ส่วนอันนี้ก็เป็น code ที่เขียนเพื่อเรียกใช้ function main ของโปรแกรม LongLexTo ที่คอมไพล์แล้วครับ library(rJava) .jinit() .jaddClassPath(“D:\\works\\LongLexTo\\R”) # replace with the path to the folder containing

ผมถูกถามว่า strain ขอเชื้อมาลาเรียจากข้อมูลที่ใช้ในโมเดลที่ทำอยู่มีอิทธิพลอย่างไรกับการเพิ่มขึ้นหรือลดลงของเชื้อที่นับได้ในคนไข้  คำถามนี้ทำให้คิดอยู่นานทีเดียวว่าจะตอบหรือแสดงกราฟอย่างไรดี เนื่องจากว่ากราฟของพารามิเตอร์ตัวหนึ่งในโมเดลแสดงการกระจายตัวเป็นแบบ bimodal ซึ่งมันอาจเกี่ยวข้องกับ strain ของเชื้อหรือเปล่าที่ทำให้มีการกระจายอย่างนั้น เพื่อที่จะตอบคำถามนี้ผมได้plot กราฟที่เรียกว่า bee swarm ครับ มันเป็นกราฟที่คล้ายๆ กับ scatter plot กับ density plot รวมกัน โดยที่เราสามารถเห็นได้ว่าจุดข้อมูลของเรานั้นมีส่วนช่วยในการสร้างกราฟอย่างไร มีการกระจายตัวอย่างไร ตามนี้ครับ ซึ่งจากกราฟผมสามารถเห็นได้ว่าข้อมูล strain จากเชื้อมาลาเรียชนิดต่างๆไม่ได้กระจุกตัวอยู่ที่ได้ที่หนึ่งนั่นก็หมายความว่าโมเดลทำนายว่า strain ไม่ได้มีอิทธิพลเกี่ยวข้องกับการการเพิ่มขึ้นหรือลดลงของเชื้อที่นับได้ในคนไข้กลุ่มนี้ กราฟลักษณะนี้สามารถ plot ได้ในโปรแกรม R โดยอาศัย library ที่ชื่อ beeswarm ครับ หรือถ้าอยากจะทำใน Mathematica สามารถดูได้ที่นี่ครับ http://mathematica.stackexchange.com/questions/42585/implementing-a-beeswarm-plot-in-mathematica