maemod package

จากที่เคยเขียนไว้เกี่ยวกับ 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]
 # d/dt (A[n]) = k[n-1]*A[n-1]
 #
 # init A[1] = 100
 # init A[2..n]=0
 #
 # k[1..n] =1
 # n = 14
 #
 #

arrayEx <-'
 !ExtraFunctions
 n<-14

!Equations
 dA <-rep(0,n)
 dA[1] <- -k[1]*A[1]
 #ArrayB
 "dA[i.indx] <- k[i.indx-1] *
 A[i.indx-1]-k[i.indx]*A[i.indx]"%@@%list("i.indx"%=>%2:(n-1))
 #ArrayE
 dA[n] <- k[n-1]*A[n-1]

!Parameters
 1,1,1,1,1,1,1,1,1,1,1

!Inits
 100,0,0,0,0,0,0,0,0,0,0,0,0,0

!Outputs
 c(dA)

!Plots

!MAEMOD_End
 '
 out<-maemod.ode(input.text = arrayEx ,timegrid = seq(0,20,0.02),sys.template = Maemod_Array)
 plot(out[,c(1,2)],type = 'l')
 cols<-rainbow(14)
 for(i in 3:n) lines(out[,c(1,i)],col=cols[i])

ในเบื้องต้นนี้ ตัวstate กับ parameter จะต้องเป็น A กับ k เท่านั้น จากตัวอย่างเราสามารถที่จะdefine ตัวแปรที่ index ต่างได้ เช่น ต้องการ

d/dt (A[2..n-1]) = k[i-1]*A[i-1]-k[i]*A[i]

โดยที่ i มีค่าตั้งแต่ 2 ถึง  n

เราก็สามารถพิมพ์ตามนี้ได้

#ArrayB

"dA[i] <- k[i-1]*A[i-1]-k[i]*A[i] "%@@%list('i'%=>%2:n)

#ArrayE

เจ้าเครื่องหมาย %@@% เป็นตัวบอกว่าเราจะใช้ array กับ สมการด้านซ้ายมือของเคื่องหมายและindexจะอยู่ด้านขวาของเครื่องหมาย ซึ่งจะต้องเป็นตัวแปรแบบ list

ส่วนเครื่องหมาย %=>% บอกว่าสัญลักษณ์หรือรูปแแบบindexด้านขวาจะมีค่าเปลี่ยนไปตาม vector ด้านซ้าย เช่น ‘i’%=>%2:n จะหมายถึง i มีค่าตั้งแต่ 2 ถึง n

ส่วน #ArrayB และ #ArrayE เป็น block ที่บอกว่ามีการใช้ ตัวแปรแบบ array ในข้อความที่อยู่ใน block นี้ครับ

codeที่มีarray blockนี้เวลาเอาไปใช้กับ maemod.ode จะต้องใช้กับtemplate เฉพาะที่ชื่อ Maemod_Array นั่นก็คือต้องเซ็ต option ชื่อ sys.template = Maemod_Array ครับ

สนใจอยากลองใช้งานก็ดูได้ที่

https://github.com/slphyx/maemod

 

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/ เรื่องที่เราจะหาค่า \alpha, \beta จากข้อมูลที่เรามีโดยที่โมเดลของเราเป็นสมการเส้นตรง y = \alpha+\beta x

อันนี้ตัวอย่าง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, 1:100)
plot(ob)

# using lm to fit the fake data
lmfit<-lm(ob$y ~ ob$time)
cat(paste0("alpha = ",format(round(lmfit$coefficients[1],4),nsmall = 4), " beta = ", format(round(lmfit$coefficients[2],4),nsmall = 4)))
abline(lmfit, col='red')

ซึ่งจากข้อมูลที่สร้างมานี้เรามารถใช้คำสั่ง lm หาได้ว่าค่า \alpha, \beta เป็นเท่าไหร่โดยคำตอบที่ได้นี้มาจากเทคนิคที่เรียกว่า QR decomposition ซึ่งจะต่างจากเทคนิค Markov chain Monte Carlo  ที่ใช้ใน Stan ที่คำตอบจะออกมาเป็นในลักษณะของdistributionที่เราประมาณได้ว่าค่าที่เป็นไปได้อาจอยู่แถวค่า median ของdistributionที่ได้นั้น

อันนี้เป็นตัวอย่างcode ที่ผมใช้fitกับข้อมูลด้านบนครับ

library(rstan)
#เป็นfunctionที่สร้างสำหรับให้ค่าเริ่มต้นของตัวแปร
initf<-function(){
 list(alpha=0.1,beta=0.1,sig=0.1)
}
#สร้างตัวแปรdataใหม่ในรูปแบบของlist
dat <- list(y=ob$y,x=ob$time,np=length(ob$y))

#codeของstan
mymodel <- "
functions{
 real fn(real x, real alpha, real beta){
 return alpha + beta*x;
 }
}
data{
 int np;
 vector[np] y;
 vector[np] x;
}
parameters{
 real alpha; 
 real beta;
 real sig;
}
model{
 vector[np] md;

alpha ~ uniform(-5,5);
 beta ~ uniform(-5,5);
 sig ~ uniform(0,30);

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

y ~ normal(md,sig);
}
generated quantities{
 vector[np] y_pred;
 vector[np] md;

for(i in 1:np)
 md[i] <- fn(x[i],alpha,beta);

for(i in 1:np)
 y_pred[i] <- normal_rng(md[i], sig);

}
"
stanout <- stan(model_code = mymodel, data = dat, chains = 1, iter = 10000, init = initf)
summary(stanout,pars=c("alpha","beta"),probs=c(0.025,0.5,0.975))$summary

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

บางคนพอเห็นcodeของ Stanแล้วบอกว่าจะมาเสียเวลาเขียนทำไมยาวๆใช้ lm สิบันทัดเดียวจบเลย ก็จริงครับสำหรับปัญหาง่ายๆแบบนี้ แต่ถ้าปัญหาที่มันซับซ้อนกว่านี้ล่ะ ผมหมายถึงโมเดลที่มันมีหลายตัวแปรที่ซับซ้อนในการคิดและแน่นอนไม่ได้เป็นแบบ linear ล่ะ Stan มันถูกออกแบบมาให้ใช้กับปัญหาที่หลากหลายครับ แต่มันก็ต้องเปลี่ยนวิธีมองปัญหาให้เป็นแบบเบย์ (Bayesian) ครับ อย่างเช่นในปัญหานี้ ….

เดี๋ยวมาต่อ

 

ใช้ประโยชน์จาก 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

##lorenz.stan
functions{
  vector lorenz(real t,vector y,vector parms){
    vector[3] ydot;
    // parms[1]=beta ,[2]=sigma,[3]=rho
    
    ydot[1]=parms[2]*(y[2]-y[1]);
    ydot[2]=y[1]*(parms[3]-y[3])-y[2];
    ydot[3]=y[1]*y[2]-parms[1]*y[3];

    return ydot;
  }
}
model{
}

จากนั้นก็ใช้ expose_stan_functions เรียกจากใน R ได้เลยตามนี้ครับ

library(rstan)

expose_stan_functions("lorenz.stan")

# parms[1]=beta ,[2]=sigma,[3]=rho
parameters1<-c(2.66666,10,28)

# x=-10, y=-12, z=30.05
init.state<-c(-10,-12,30.05)

nderivs<-function(t,state,parameters){
  tmp<-lorenz(t,state,parameters)
  return(list(tmp))
}

run1<-ode(y=init.state,times=seq(0,15,0.01), 
   func=nderivs, parms=parameters1)

plot(run1)

 

 

word cloud ภาษาไทยใน R

ตัวอย่างการทำเมฆกลุ่มคำใน 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]])

พอได้vectorของแต่ล่ะคำแล้วก็สร้าง data frame สำหรับนับว่าแต่ล่ะคำนั้นมีความถี่หรือมีจำนวนเท่าใด

wordfreq<-as.data.frame(table(outtxt))

จากนั้นก็เอาตัวแปรสำหรับเก็บความถึ่ของคำนี้ไปสร้างwordcloudได้เลยครับ

wordcloud2(wordfreq,color = "random-light",shape = 'star',backgroundColor = "black")

star

ส่วนอันนี้ก็เป็นอีกตัวอย่างครับ โดยที่ไฟล์ข้อความผมcopyมาจาก http://king.kapook.com/royal_words_2548.php

RLongLexToF(inputfilename = "test2.txt",outputfilename = "outtest2.txt")
outtxt2<-as.vector(strsplit(readLines("outtest2.txt"),"[|]")[[1]])

wordfreq2<-as.data.frame(table(outtxt2))
wordcloud2(wordfreq2)

cloud2

ความยาวของตัวอักษรที่ R console จะรับได้ในหนึ่งคำสั่ง

กำลังเขียน 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