Category: R

มิกกี้เม้าส์ใน R

วาด Mickey mouse ใน R

library(showtext)
library(ggplot2)

font_add("TH Mali Grade 6", 
"PATH_to_/TH Mali Grade6.ttf")

showtext_auto()
# Define the functions for x(t) and y(t)
x_func <- function(t) {
  -119/9 * sin(1/10 - 25 * t) - 945/43 * sin(5/11 - 23 * t) - 501/8 * sin(1/47 - 10 * t) - 
    801/20 * sin(13/10 - 8 * t) - 607/10 * sin(7/10 - 7 * t) - 1975/9 * sin(3/7 - 2 * t) + 
    4655/8 * sin(t + 29/10) + 773/7 * sin(3 * t + 76/17) + 1138/9 * sin(4 * t + 1/3) + 
    619/4 * sin(5 * t + 5/8) + 3355/26 * sin(6 * t + 13/10) + 595/9 * sin(9 * t + 28/11) + 
    905/9 * sin(11 * t + 40/11) + 149/10 * sin(12 * t + 7/11) + 153/5 * sin(13 * t + 2) + 
    127/8 * sin(14 * t + 13/6) + 225/8 * sin(15 * t + 3/2) + 49/4 * sin(16 * t + 25/7) + 
    133/9 * sin(17 * t + 15/8) + 301/15 * sin(18 * t + 31/9) + 232/7 * sin(19 * t + 11/7) + 
    438/19 * sin(20 * t + 10/7) + 67/3 * sin(21 * t + 19/10) + 114/11 * sin(22 * t + 38/37) + 
    153/11 * sin(24 * t + 10/9) + 153/8 * sin(26 * t + 11/12) + 443/34 * sin(27 * t + 23/14) + 
    39/7 * sin(28 * t + 45/11) + 32/3 * sin(29 * t + 4/11) + 110/13 * sin(30 * t + 3/8) + 
    41/7 * sin(31 * t + 70/23) + 39/7 * sin(32 * t + 13/9)
}

y_func <- function(t) {
  -86/19 * sin(6/7 - 30 * t) - 89/8 * sin(2/3 - 29 * t) - 687/49 * sin(4/9 - 24 * t) - 
    294/11 * sin(11/8 - 20 * t) - 135/4 * sin(10/11 - 17 * t) - 71/7 * sin(1/6 - 15 * t) - 
    226/5 * sin(3/8 - 14 * t) - 295/8 * sin(1/8 - 11 * t) - 220/7 * sin(17/16 - 10 * t) - 
    708/7 * sin(1/6 - 9 * t) - 2017/8 * sin(5/4 - 2 * t) - 1240/3 * sin(11/10 - t) + 
    104/15 * sin(18 * t) + 2216/7 * sin(3 * t + 39/11) + 601/3 * sin(4 * t + 18/11) + 
    468/11 * sin(5 * t + 34/11) + 127/7 * sin(6 * t + 27/7) + 457/8 * sin(7 * t + 25/6) + 
    201/4 * sin(8 * t + 23/7) + 71/3 * sin(12 * t + 4/7) + 340/7 * sin(13 * t + 39/11) + 
    216/7 * sin(16 * t + 27/13) + 31/5 * sin(19 * t + 24/25) + 8/5 * sin(21 * t + 1/7) + 
    217/11 * sin(22 * t + 23/14) + 34/3 * sin(23 * t + 19/8) + 157/9 * sin(25 * t + 3/10) + 
    69/4 * sin(26 * t + 2/3) + 29/5 * sin(27 * t + 15/7) + 58/7 * sin(28 * t + 29/8) + 
    95/11 * sin(31 * t + 19/8) + 34/9 * sin(32 * t + 40/9)
}

# Generate values
t_values <- seq(0, 2 * pi, length.out = 1000)
x_values <- sapply(t_values, x_func)
y_values <- sapply(t_values, y_func)

df <- data.frame(t = t_values, x = x_values, y = y_values)

# Plot using ggplot2
p <- ggplot(df, aes(x = x, y = y)) + 
  geom_path(color = "black", size = 1) + 
  theme_minimal() + 
  geom_text(aes(x = 500, y = 800, label = "กลุ่มผู้ใช้ R"), 
            family = "TH Mali Grade 6", color = "black", size = 26)

print(p)

เรียน R Shiny

สำหรับผู้ที่สนใจจะพัฒนาเว็บและแอปใน R ด้วย Shiny ครับ เนื้อหาที่ผมสอนใน workshop ซึ่งเป็นภาษาไทย ได้ถูกเผยแพร่ไว้ที่นี่ Shiny เบื้องต้น (ammnet-thailand.github.io) ครับ

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

%d bloggers like this: