a simple pendulum in R

เป็นตัวอย่างที่ผมเขียนขึ้นตอนใช้สอนการเขียน function ในภาษา R ครับ

library(magick)
# the formula for the changing of angle over time 
# https://en.wikipedia.org/wiki/Pendulum
theta <- function(t,period,theta0){
  theta0*cos(2*pi*t/period)
}

# for drawing the pendulum at a given point (x,y), i.e. drawing a line between (0,0) and (x,y)
draw.pendulum <- function(pen.x, pen.y){
  
  pivot.x <- 0
  pivot.y <- 0

  plot(x=pivot.x,y=pivot.y,type = 'l',ylim = c(-1.5,0.1),asp = 1)
  lines(x=c(pivot.x,pen.x),y=c(pivot.y,pen.y))
  points(x = pen.x, y=pen.y,col='red')
}

# calculate the position of the pendulum over time
# l-length of the pendulum
# theta0 - the starting angle (in radian)
# runtime - the simulation time (in seconds)
# step - time step
pendulum <- function(theta0,l,runtime,step = 0.1){
  g <- 9.8
  period <- 2*pi*sqrt(l/g)
  
  times <- seq(0.0,runtime, by = step)

  pen.x <- l*sin(theta(times, period = period,theta0 = theta0)) 
  pen.y <- -l*cos(theta(times,period = period,theta0 = theta0))
  
  output <- matrix(c(times,pen.x,pen.y),ncol = 3)
  return(output)
}


pendulum.sim <- function(theta0, l, runtime, step = 0.1){
  
  sim.data <- pendulum(theta0 = theta0, l=l, runtime = runtime)
  
  ## export each frame as png
  for(i in 1:nrow(sim.data)){
    png(filename = paste0(tempdir(),"/",i,".png"))
    draw.pendulum(sim.data[i,2],sim.data[i,3])
    dev.off()
  }
  
  frames <- image_read(paste0(tempdir(),"/",1:nrow(sim.data),".png"))
  animate <- image_animate(image = frames, fps = 1/step)
  
  return(animate)
  
}

animate <- pendulum.sim(theta0 = pi/4, l = 1, runtime = 10)
print(animate)

มาลองใช้ Stan แก้ปัญหา

พอดีว่าง ระหว่างรอผลคำนวณอะไรบางอย่าง ผมเห็นมีคนโพสท์ถามที่เวบพันทิปตามนี้ครับ

จากคำถามนี้เราสามารถที่จะใช้ Stan แก้ปัญหาได้ถ้าผมมองว่าการวัดเปรียบเทียบเครื่องมือมาตรฐาน(จริง)กับที่ปรับปรุงขึ้นมานั้นในแต่ล่ะครั้งนั้นไม่ได้เกี่ยวกันเลย และค่าจากวัดของเครื่องที่ปรับปรุงนั้นมีการกระจายรอบค่าจากเครื่องมาตราฐานแบบ normal distribution โดยมีค่า standard deviation หรือ SD อยู่ค่าหนึ่งที่เป็นตัวกำหนดประสิทธิภาพของเครื่อง อย่างเช่น จากข้อมูลที่ให้มาเมื่อวัดเทียบกับเครื่องมือจริงที่วัดได้ 2 แต่ค่าจากเครื่องปรับปรุงวัดมา 3 ครั้งได้ (1.8,2.1,1.9) เราจะสมมุติให้ทั้ง 3 ค่านี้กระจายรอบค่าใดค่าหนึ่ง โดยที่ ถ้าเราปรับปรุงได้เจ๋ง ค่าที่มันกระจายรอบนี้มันก็ควรจะได้เท่ากับค่าจากเครื่องจริงโดยที่มีค่า SD น้อยที่สุดเท่าที่จะเป็นไปได้ แต่ในความเป็นจริงก็จะมีผลผันผวนหรือerrorต่างๆเข้ามาเกี่ยวข้อง ที่ผมจะทำให้ดูนี้เราจะลองใช้โปรแกรมอย่าง Stan มาช่วยหาดูว่าค่ากลางที่เครื่องปรับปรุงหรือสร้างขึ้นมานั้นมันกระจายรอบในแต่ล่ะครั้งของการวัดมันคือค่าอะไรและมี SD เท่าใด

ผมใช้โปรแกรมที่ชื่อ Stan ในภาษา R ผ่าน library ที่เรียกว่า rstan ครับ ด้านล่างนี้ก็เป็น code ที่ผมใช้กับปัญหานี้ครับ

#โหลด library ที่จะใช้งานครับ 
#ในที่นี้มี 2 ตัวคือ rstan กับ bayesplot (ใช้วาดกราฟสรุปผล)
library(rstan)
library(bayesplot)

#เตรียมข้อมูล โดยผมแยกตามค่าของเครื่องจริง
data <- list(
  ob2 = c(1.8,2.1,1.9),
  ob3 = c(3.2,3.0,3.2),
  ob4 = c(3.9,4.3,4.4),
  ob5 = c(5.0, 5.2, 5.5),
  ob6 = c(6.1,5.9,6.5),
  n=3
)

# code ของโมเดล
model <- "
data{
      //จำแหนกประเภทของตัวแปรของข้อมูล ว่าเป็นเลขจำนวนเต็มหรือทศนิยม พร้อมกำหนดขนาด
        int<lower=1> n; 
	real ob2[n];
	real ob3[n];
	real ob4[n];
	real ob5[n];
	real ob6[n];
}
parameters{
// กำหนดประเภทของตัวแปรที่จะใช้ในโมเดล ซึ่งในที่นี้คือค่ากลางที่ข้อมูลมันกระจายรอบ ซึ่งแบ่งตามค่าจริง
// sig เป็นค่า SD 
	real ob2mu;
	real ob3mu;
	real ob4mu;
	real ob5mu;
	real ob6mu;
	real<lower=0> sig;
}
model{
//กำหนดว่าค่ากลางมันและSD อยู่ในช่วงไหน จากการกระจายแบบไหน 
//ซึ่งในที่นี้ผมให้มันมาจาก uniform distribution โดยใส่ช่วงที่คิดว่าค่ามันจะอยู่ในนั้น
	sig ~ uniform(0,3);
	ob2mu ~ uniform(0,7);
	ob3mu ~ uniform(0,7);
	ob4mu ~ uniform(0,7);
	ob5mu ~ uniform(0,7);
	ob6mu ~ uniform(0,7);

// กำหนดให้ค่าที่วัดมาในแต่ล่ะค่าจริงกระจายแบบ normal รอบค่ากลางอันหนึ่ง
	ob2 ~ normal(ob2mu, sig);
	ob3 ~ normal(ob3mu, sig);
	ob4 ~ normal(ob4mu, sig);
	ob5 ~ normal(ob5mu, sig);
	ob6 ~ normal(ob6mu, sig);
 
}
"
fit <- stan(model_code=model, data = data, iter = 10000)

ผลลัพท์ที่ได้ก็จะประมาณนี้ครับ

เห็นได้ว่าค่ากลางหรือ median ที่ 50% นั้นค่อนข้างจะใกล้กับค่าจริงคือ เช่นที่ค่าจริง = 2 เครื่องปรับปรุงทำได้ 1.93(95%CI: (1.63,2.24)) หรือที่ 6 เครื่องปรับปรุงทำได้ที่ 6.15 (CI: (5.86,6.48)) และ SD = 0.25 CI: (0.17-0.44)

posterior <- as.matrix(fit)
plot_title <- ggtitle("Posterior distributions",
                      "with medians and 95% intervals")
mcmc_areas(posterior,
           pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ),
           prob = 0.95) + plot_title
ตัวอย่างการกระจายของค่าที่ได้
posterior2 <- extract(fit, inc_warmup = T, permuted = FALSE)

color_scheme_set("mix-blue-pink")
p <- mcmc_trace(posterior2,  pars = c("ob2mu","ob3mu","ob4mu","ob5mu","ob6mu", "sig" ), 
                n_warmup = 5000,
                facet_args = list(nrow = 6, labeller = label_parsed))
p + facet_text(size = 15)

จากกราฟการกระจายที่ได้จะเห็นได้ว่าผมใส่เป็น uniform distribution ที่ช่วงระหว่าง 0-7 เลยแต่ผลลัพท์ที่ได้เป็น normal distribution ครับ

chains ของแต่ล่ะตัวแปร

Kung package แพ็คเกจกุ้ง

แพ็คเกจกุ้งเป็นแพ็คเกจภาษา R ที่ผมเขียนขึ้นมาสำหรับใช้สร้าง Shiny App สำหรับโค้ดของโมเดลที่สร้างจากสมการอนุพันธ์หรือ ODE ที่ใช้ตัว solver จากแพ็คเกจที่ชื่อ deSolve ครับ

ไอเดียก็มีเพียงว่าจากสมการ ode ที่สร้างขึ้นมาด้วยถ้าตัวแปรในสมการมีการเปลี่ยนแปลงค่า ผลลัพท์ที่ออกมามันจะเปลี่ยนแปลงไปอย่างไร การเขียนใน R นั้นเราก็อาจจะใช้คำสั่ง manipulate มาช่วยได้แต่มันก็ยังมีข้อจำกัดเรื่องจำนวนของตัวแปรที่สามารถเรียกมาใช้ได้ เพราะมันขึ้นอยู่กับขนาดหน้าจอของตัว plot ผมเลยคิดว่าถ้าใช้ shiny มันน่าจะดูดีกว่าและสามารถใส่ตัวแปรได้มากกว่า ขึ้นกับการออกแบบ ui แต่การเขียน shiny ก็ไม่ใช่ว่าจะง่ายสำหรับคนเริ่มเรียนรู้ R ผมเลยคิดว่าน่าจะเขียนแพ็คเกจที่สามารถช่วยเรื่องนี้ได้

แพ็คเกจที่เขียนก็พยายามทำให้มันใช้ง่ายมากที่สุดจาก code ที่มีอยู่แล้ว โดยผู้ใช้เพิ่มคำสั่งไม่มากก็สามารถที่จะสร้าง shiny application ได้แล้ว ตัวอย่างเช่น ถ้าผู้ใช้มี code อยู่แล้วตามนี้ ซึ่งมันจะคำนวณและ plot กราฟให้ตามค่าของตัวแปรพารามิเตอร์ที่กำหนดไว้ เช่น gamma = 0.14286 และ beta = 0.6

library(deSolve)

init <- c(S = 1-1e-6, I = 1e-6, R= 0.0)
times <- seq(0, 70, by = 1)

sir <-function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I
    
    return(list(c(dS, dI, dR)))
  })
}

parameters <- c(
  gamma = 0.14286, 
  beta = 0.6
)

out <- as.data.frame(ode(y = init, times = times, func = sir, parms = parameters))

matplot(times, out[,c("S","I","R")], type = "l", xlab = "Time", ylab = "Susceptibles and Recovereds", main = "SIR Model", lwd = 1, lty = 1, bty = "l", col = 2:4)
legend(40, 0.7, c("Susceptibles", "Infecteds", "Recovereds"), pch = 1, col = 2:4)

ถ้าผู้ใช้ต้องการดูว่ากราฟจะเปลี่ยนแปลงอย่างไรถ้าตัวแปร gamma หรือ beta มีการเปลี่ยน ก็เพียงเพิ่ม code ลงไปในส่วนต่างๆ ตามนี้

!Start

library(deSolve)

init <- c(S = 1-1e-6, I = 1e-6, R= 0.0)
times <- seq(0, 70, by = 1)

sir <-function(time, state, parameters) {
  with(as.list(c(state, parameters)), {
    dS <- -beta*S*I
    dI <- beta*S*I - gamma*I
    dR <- gamma*I

    return(list(c(dS, dI, dR)))
  })
}

!Parameters
parameters <- c(
  gamma = 0.14286,
  beta = 0.6
)

!ODECMD
out <- as.data.frame(ode(y = init, times = times, func = sir, parms = parameters))

!PostProcess

!Plots
matplot(times, out[,c("S","I","R")], type = "l", xlab = "Time", ylab = "Susceptibles and Recovereds", main = "SIR Model", lwd = 1, lty = 1, bty = "l", col = 2:4)
legend(40, 0.7, c("Susceptibles", "Infecteds", "Recovereds"), pch = 1, col = 2:4)

!Controls
sliderInput("beta","beta", min = 0,max = 10,step = 0.001,value = 1.4),
sliderInput("gamma","gamma", min = 0,max = 1,step = 0.001,value = 0.14)

!End

จะเห็นว่ามี keywords ที่เพิ่มเข้าไปดังนี้

!Start – เป็นตัวบอกเริ่มต้น เราสามารถบอกรายละเอียดของตัวแปรอย่าง เวลา และค่าเริ่มต้นของ state หรือ compartment ต่างๆและ function ของระบบสมการหลักที่จะแก้ด้วย deSolve::ode ใน keywords นี้ได้เลยครับ

!Parameters – เป็นตัวแปร parameters ของโมเดล ode ครับ ซึ่งต้องกำหนดเป็น vector ตามตัวอย่างนี้เลยครับ

!ODECMD – สำหรับคำสั่งที่ต้องใช้คำสั่ง ode จาก deSolve ครับ สำหรับตอนนี้ต้องให้ตัวแปร output ชื่อ out นะครับ

!PostProcess – ถ้ามีการคำนวณค่าอะไรบางอย่างจาก out ก็ให้ใส่ใน keyword นี้ได้เลยครับ

!Controls – สำหรับใส่ sliders ของตัวแปรที่สนใจครับ

!End – สำหรับบอกว่าจบแล้ว 🙂

โดยหลังจากที่ใส่ keywords และsave เป็นไฟล์ใหม่แล้ว เช่นชื่อ mysystem.R เราสามารถสร้าง shiny app ได้ด้วยการรันคำสั่งนี้ครับ

runSystem(‘mysystem.R’)

ผลที่ได้ก็จะประมาณนี้ครับ

ถ้าสนใจอยากใช้งานก็ทำการติดตั้งได้โดยพิมพ์

devtools::install_github('slphyx/Kung')

ดูเพิ่มเติมที่ https://github.com/slphyx/Kung

แสดง output จาก ggplot ใน Mathematica

ลองดูตัวอย่างที่ผมเขียนนี้ครับ  ไอเดียก็คือว่าให้มันเขียนภาพเป็นไฟล์ pdf ด้วยคำสั่ง ggsave แล้วเรียกกลับมาใน Mathematica ด้วย Import  โค้ดที่เขียนนี้ดัดแปลงมาจาก https://mathematica.stackexchange.com/questions/16726/how-to-get-a-plot-generated-by-r-returned-in-an-output-cell

(* using R that already installed on the computer *)
R363 = "C:\\Data-Work\\Programs\\R-3.6.3\\App\\R-Portable";
InstallR["RHomeLocation" -> R363, "RVersion" -> "3.6.3", 
 "NativeLibLocation" -> 
  "C:\\Data-Work\\Programs\\R-3.6.3\\App\\R-Portable\\library\\rJava\\\
jri\\x64"]

REvaluate["
 library(ggplot2)
 "]
(* export ggplot as pdf file *)
Wrapper = RFunction["function(filename,plotfun){
   ggsave(filename, plot=plotfun(), width=6, height=4, \
device=cairo_pdf)
   }"];

(* generate the plot as pdf and import back to show in Mathematica *)


Getggplot[plotFun_] := Module[{tempfile, rcode, runcode},
   tempfile = FileNameJoin[{$TemporaryDirectory, "temp.pdf"}];
   If[FileExistsQ[tempfile], Quiet@DeleteFile[tempfile]];
   rcode = "ADA <- ggplotcmd;";
   runcode = StringReplace[rcode, {"ggplotcmd" -> plotFun}];
   Wrapper[tempfile, RFunction["function()" <> "{" <> runcode <> "}"]];
   If[! FileExistsQ[tempfile], Return[$Failed]];
   Import[tempfile] // First
   ];

RDotNet 1.7 ใช้งานได้กับ R.3.3-4 เท่านั้น

ผมเขียนโปรแกรมเล็กๆที่เรียกใช้งาน RDotNet ไว้หลายตัว ซึ่งมันก็ทำให้สะดวกมากๆในการที่จะเรียก R มาใช้งานใน C# ซึ่งก่อนหน้านี้จะใช้งานกับ R ที่ค่อนข้างเก่าหน่อยพวก R 3.3-4 แต่พอมี library บางตัวที่บังคับให้ต้อง upgrade ตัว R ไป version ที่สูงขึ้นอย่าง 3.5 นี้ทำให้โปรแกรมที่เขียนหลายตัวมันมีปัญหาทันที เหตุผลก็เพราะมัน GetInstance() ไม่ได้ มีการพูดถึงปัญหานี้ที่ https://github.com/jmp75/rdotnet/issues/70 ตอนนี้ที่ทำได้ก็คือถ้า library ที่อยากใช้นั้นถ้าไม่ได้ใช้ features ที่มันมีมาใหม่นั้นจริงๆก็ต้องใช้ R ตัวเก่าไปก่อน

 

%d bloggers like this: