เรื่องปวดหัว 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

 

การ 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

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

LOWESS ใน Mathematica

ตั้งแต่มี RLink มานี่การทำอะไรที่เกี่ยวกับทางสถิติใน Mathematica มันดูง่ายไปหมด แต่ก่อนเวลาผมทำพวก local regression ผมต้อง export ข้อมูลจาก Mathematica ไปทำใน R หรือ Stata แต่เดี่ยวนี้ไม่ต้องแล้ว ก็ทำมัน ใน Mathematica นี่แหละด้วยการใช้ RLink เรียกคำสั่ง lowess (stats) จาก R อีกที

อันนี่เป็นตัวอย่างที่ผมใช้ทำ LOWESS (http://en.wikipedia.org/wiki/Local_regression)ครับ ผมมีข้อมูลที่เป็นคู่อันดับ (x,y) โดย x ของผมก็คือเวลา และ y คือข้อมูลจากการวัด วิธีการที่ทำก็คือสร้าง function ที่ชื่อ Lowess โดยให้มันอ่านค่า List ซึ่งเป็นคู่อันดับของข้อมูลแล้วแปลงเป็น Data Frame เพื่อเอาไปใช้กับคำสั่ง lowess ใน R ครับ

lowess

ตัวอย่างที่ได้ครับ

gam1

ตัวอย่างการใช้งาน RLink ใน Mathematica

ตัวอย่างการใช้โปรแกรมทางสถิติ R ใน Mathematica ผ่าน RLink ครับ

**รายละเอียดต่าง ๆ ของการใช้งานของ RLink สามารถดูเพิ่มเติมได้จาก Document center ของ Mathematica ครับ

RLink เป็นแพ็คเกจอันหนึ่งของ Mathematica (สำหรับ version 9 ขึ้นไปครับ) ที่ช่วยให้เราสามารถที่จะ Run คำสั่งของ R หรือโปรแกมที่เราเขียนขึ้นด้วย R  ได้โดยตรงจากในตัวของ Mathematica ครับ อีกทั้งเรายังสามารถที่จะส่งผ่านหรือรับข้อมูล/ผลลัพท์ระหว่าง Mathematica กับ R ได้โดยตรง

R เป็นโปรแกรมทางสถิติที่ได้รับความนิยมสูงมากทั้งในการเรียนการสอนและการวิจัย มีคำสั่งหรือชุดคำสั่ง (library/package) สำหรับทำงานเฉพาะด้านที่เราสามารถนำมาประยุกต์ใช้ในงานของเราได้ (http://cran.r-project.org/web/packages/available_packages_by_name.html)

การใช้งาน R ผ่าน RLink ใน Mathematica นี้จะเป็นการเพิ่มขีดความสามารถของ Mathematica ในด้านสถิติให้มากขึ้นครับ

ตัวอย่างการใช้งาน RLink ครับ

ก่อนอื่นเลยเราต้องทำการโหลด RLink เข้ามาใน Mathematica ก่อนครับ

<<RLink`

จากนั้นก็ติดตั้งตัวโปรแกรม R ที่ทาง Wolfram เตรียมไว้ให้แล้วด้วยคำสั่ง

RLinkResourcesInstall[]

การติดตั้งตัวโปรแกรม R ด้วย RLinkResourcesInstall นี้ทำครั้งแรกครั้งเดียวพอครับ การใช้งานครั้งต่อๆไปไม่ต้องครับ

เวลาเรียกใช้งาน R ก็พิมพ์

InstallR[]

แต่ถ้าใครอยากใช้กับ R ที่ติดตั้งเรียบร้อยแล้วในเครื่องของตัวเอง (เช่นที่ C:\\Program Files\\R\\R-2.15.1) ก็เพียงพิมพ์

InstallR[“RHomeLocation” -> “C:\\Program Files\\R\\R-2.15.1”]

คำสั่งหลัก ๆ ใน RLink จะพูดถึงนี้ก็มี REvaulate, RSet และ RFunction

REvaluate ดูเหมือนจะเป็นพระเอกของ RLink ครับ มันใช้สำหรับการรันคำสั่งต่าง ๆ ของ R พร้อมกับรับค่าที่ได้ส่งกลับมาที่ Mathematica ครับ (กรณีที่คำสั่ง R นั้นไม่มี ; ปิดท้าย) เช่น

Revaluate[“{

a <- 2;

b <- 1.334;

(a*a)-b

}”]

RSet เป็นคำสั่งที่ช่วยให้เราส่งค่าพวกตัวเลขหรือลิสต์ จาก Mathematica ให้กับตัวแปรใน R

เช่นต้องการให้ตัวแปร x ที่เราใช้ใน R มีค่าเท่ากับลิสต์ของเลขสุ่มระหว่าง 0 ถึง 1 ที่เรียงจากมากไปน้อย

RSet[“x”, Sort[#, Greater] &@RandomReal[1,10]]

 

RFunction  เป็นคำสั่งที่ช่วยให้เราเรียกใช้ function ที่เราเขียนใน R มาใช้ใน Mathematica ซึ่งวิธีการเรียกก็เหมือนกับ function หรือคำสั่งของ Mathematica เลยครับ เช่น  ถ้าผมต้องการ function ชื่อ stats ที่คำนวณค่า mean,median,sd และ variance ของlist ของตัวเลข ผมก็สามารถเขียนได้ตามนี้เลยครับ

stats = RFunction[“function(x){ c(mean(x),median(x),sd(x),var(x)) } “]

stats[RandomReal[1,100]]

 

ตัวอย่างที่ผมลองเอาไปใช้งานครับ ผมมีข้อมูลที่เป็นคู่อันดับอยู่หลายชุด และแต่ล่ะชุดจะมีคู่อันดับหรือจุดไม่เกิน 7 จุดและผมได้ทดลอง fit (ทำ linear regression) กับ เส้นตรง y = ax+c ไปแล้ว แต่พอดูด้วยตาแล้วสมการกำลังสอง y = ax^2+bx+c น่าจะ fit ได้ดีกว่าเส้นตรงในบางชุดของข้อมูล ก็เลยคิดว่าน่าเปรียบเทียบทั้งสองสมการนี้ด้วย Likelihood ratio test โดยใช้ p-value เป็นตัวตัดสินใจครับ ผมเลือกเอาว่าถ้า p-value มันน้อยกว่า 0.05 แสดงว่าสมการกำลังสองมัน fit ได้ดีกว่าเส้นตรง และถ้ามันมากกว่าหรือเท่ากับ 0.05 แสดงว่าสมการเส้นตรงมัน fit ได้ดีกว่า  ผมใช้คำสั่ง lrtest ของpackageชื่อ lmtest ในการทำ likelihood ratio test ครับ

แต่จาก R ที่ติดตั้ง ด้วยคำสั่ง RLinkResourcesInstall[] จะไม่มีแพ็คเกจ lmtest มาด้วย ดังนั้นต้องติดตั้งก่อนครับ

เราสามารถที่จะติดตั้ง R package ตัวที่ต้องการเพิ่มได้โดยเพียงพิมพ์คำสั่งนี้ใน Mathematica ได้เลย

REvaluate[“install.packages(\” ชื่อ package ที่ต้องการ “\)”]

มันจะมี windows อันหนึ่ง popup ขึ้นมาให้เราเลือก mirror host สำหรับ download ไฟล์ package ครับ

ดังนั้นในกรณีของผมนี้ก็เพียงพิมพ์

REvaluate[“install.packages(\”lmtest”\)”]

ตัวอย่าง code ครับ

rlinkcode1

rlinkgraphs

ติดตั้ง R package เพิ่มใน RLink

ใครที่ใช้ R ผ่าน RLink ของ Mathematica 9.0 บน Windows สามารถที่จะติดตั้ง R package ตัวที่ต้องการได้โดยเพียงพิมพ์

REvaluate[“install.packages(\” ชื่อ package ที่ต้องการ “\)”]

มันจะมี windows อันหนึ่ง popup ขึ้นมาให้เราเลือก mirror host สำหรับ download ไฟล์ package

เวลาเรียกใช้ก็เหมือนใน R เลยครับ เช่น

REvaluate[“{

library(ชื่อ library)

blah

blah ..

}”]

 

🙂