Probabilistic Programming Language

รู้สึกว่าตั้งแต่เรื่อง Machine Learning กำลังเป็นที่พูดถึงเรื่อยๆ เรื่องของ Probabilistic Programming Language ก็เป็นที่พูดถึงมากขึ้นไปด้วย แต่บทความที่พูดถึงว่า Probabilistic Programming Language (PPL)  คืออะไร มีประโยชน์อย่างไรและเอาไปใช้งานยังไงในบ้านเรานั้นน้อยมาก ผมก็เลยคิดว่าผมน่าจะแชร์อะไรบ้างตามที่ได้ไปรู้ไปเห็นมาครับ

PPL คืออะไร PPL ก็เป็นรูปแบบของภาษาหรือชุดคำสั่งในคอมพิวเตอร์แบบหนึ่งที่จะช่วยให้เราทำสถิติเชิงอนุมาน Inference กับโมเดลความน่าจะเป็นต่างๆ ได้ง่ายขึ้น ซึ่งแน่นอนมันก็จะเกี่ยวข้องกับทฤษฎีของเบย์ด้วย  ลองดูตัวอย่างของการทำ Linear regression นี้ดูนะครับ สมมุติว่าเรามีข้อมูลเป็นคู่ (x_1,y_1),(x_2,y_2),...,(x_n,y_n) และมีโมเดลที่ใช้อธิบายข้อมูลเป็นสมการเส้นตรง y_i=\alpha+\beta x_i+\epsilon_i โดยที่ \epsilon_i \sim N(0,\sigma^2) (หรือในบางตำราจะใช้ precision หรือความถูกต้อง \tau = 1/\sigma^2) หรือ \epsilon_i = y_i-\alpha-\beta x_i เราสามารถที่จะเขียน probability density function (PDF) สำหรับ likelihood ที่เราจะได้ y_i  ได้ว่า

f(y_i|\alpha,\beta,\tau,x_i) = \frac{\sqrt{\tau}}{\sqrt{2 \pi}} \exp(-\frac{\tau}{2} (y_i-\alpha-\beta x_i)^2)

ดังนั้น likelihood สำหรับทุกๆจุดของข้อมูลก็สามารถเขียนได้ว่า (เอาของทุกจุดคูณกัน)

L(\alpha,\beta,\tau|x,y) =\prod_{i=1}^{n}f(y_i|\alpha,\beta,\tau,x_i)

= \tau^{\frac{n}{2}}(2\pi)^{-\frac{n}{2}}\exp(-\frac{\tau}{2} \sum_{i=1}^{n}(y_i-\alpha-\beta x_i)^2 )

พิจารณาที่เทอมยกกำลังสองนะครับ

\sum_{i=1}^{n}(y_i-\alpha-\beta x_I)^2 = \sum_{i=1}^{n}(y_i^2+n\alpha^2+\beta^2\sum_{i=1}^{n}x_i^2 -2\alpha n \bar{y} -2\beta\sum_{i=1}^{n}x_i y_i +2\alpha\beta n\bar{x})

โดยที่ \bar{x}=\frac{1}{n}\sum_{i=1}^{n}x_i และ \bar{y}=\frac{1}{n}\sum_{i=1}^{n}y_i

ถ้าพิจารณาดูที่ \alpha เราสามารถเขียนฟังก์ชั่นของความน่าจะเป็นของ \alpha ได้ว่า f(\alpha|\beta,\tau,x,y) \propto \exp(-\frac{n\tau}{2}(\alpha^2-2\alpha(\bar{y}-\beta\bar{x}))) หรือ \alpha นั้นสามารถสุ่มมาจาก Normal distribution ตามนี้เลย

\alpha |\beta,x,y \sim N(\bar{y}-\beta\bar{x},\frac{1}{n\tau})

มันมาจากการโยกย้ายฟ้าดิน 🙂 บวกกับอาศัยความสัมพันธ์นี้ครับ ถ้า Y มี PDF g(y) โดย g(y) \propto \exp(-\frac{1}{2}(Ay^2-2By)) เราจะได้ว่า  Y \sim N(\frac{B}{A},\frac{1}{A}) หรือ ถ้า g(y) \propto y^A \exp(-By) เราจะได้ว่า Y \sim Gamma(A+1,B) (พยายามนึกอยู่ว่าไปจำมาจากไหน)

ทำในลักษณะคล้ายๆกันกับ \beta และ \tau เราจะได้ว่า

\beta|\alpha,\tau,x,y \sim N((\sum_{i=1}^{n} x_i y_i-\alpha n\bar{x})/ \sum_{i=1}^{n}x_i^2,(\tau\sum_{i=1}^{n}x_i^2)^{-1})

\tau|\alpha,\beta,x,y \sim Gamma(\frac{n}{2}+1, \frac{1}{2}\sum_{i=1}^{n}(y_i-\alpha-\beta x_i)^2)

นั่นก็หมายความว่าในการที่เราจะหาค่า \alpha , \beta , \tau เราก็เพียงสุ่มจากสูตรที่ได้นี้ได้เลยครับ ซึ่งเราสามารถที่จะใช้ R ช่วยได้เลย ครับ ****(จริงๆแล้วในตัวอย่างนี้ ผมสมมุติให้เราไม่รู้อะไรเลยสำหรับ prior distribution ( i.e.,g(\alpha)) หรือข้อมูลอะไรบางอย่าง เช่นของ \alpha , \beta , \tau   มีรูปแบบการกระจายตัวหรือ distribution เป็นอย่างไร มีช่วงของค่าที่เป็นไปได้อย่างไร  ง่ายที่สุดผมให้  prior distribution g(\alpha)  \propto 1 หรือ \int_{-\infty}^{\infty} g(\alpha) \mathrm{d} \alpha = \infty เหมือนกันสำหรับ \beta, \tau) ครับ

\alpha |\beta,x,y \sim N(\bar{y}-\beta\bar{x},\frac{1}{n\tau}) เราเขียนใน R ได้ว่า

alpha<-rnorm(1,mean(y)-beta*mean(x),1/(n*tau))

\beta|\alpha,\tau,x,y \sim N((\sum_{i=1}^{n} x_i y_i-\alpha n\bar{x})/ \sum_{i=1}^{n}x_i^2,(\tau\sum_{i=1}^{n}x_i^2)^{-1})

m<-(sum(x+y)-alpha*n*mean(x))/sum(x**2)

s<-1/(tau*sum(x**2))

beta<-rnorm(1,m,s)

\tau|\alpha,\beta,x,y \sim Gamma(\frac{n}{2}+1, \frac{1}{2}\sum_{i=1}^{n}(y_i-\alpha-\beta x_i)^2)

w<- y-alpha-beta*x

tau<-rgamma(1,((n/2)+1),(sum(w**2)/2))

ทดลองเอาไปใช้จริงอาศัย R code ที่ผมเขียนข้างล่างนี้ดูครับ  วิธีการเอาไปใช้ก็ลองสร้างข้อมูลมาสักหนึ่งชุด ในที่นี้ผมมีfunction ชื่อfakedata สำหรับสร้างข้อมูลครับ ก็ใช้สมการเส้นตรงนี่แหละแล้วก็บวกค่า error ที่สุ่มมาบวกเข้าไป โดยในตัวอย่างนี้ผมลองสร้างโดยให้ \alpha =1.5, \beta=0.5 และก็ให้ \sigma =5 จากนั้นก็ลองใช้ คำสั่ง lm ของ R ลอง fit ดูว่าจะได้ \alpha \beta เท่าไหร่ ซึ่งค่านี่แหละเอาไว้เปรียบเทียบกับที่เราจะได้จากการสุ่มมาจากสูตรที่เพิ่งคิดไปครับ

#fake data
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')


#
# linear regression from the posterior distributions
# slphyx
LinReg<-function(y,x,size){
  #initial values
  alpha<-0
  beta<-1
  tau<-1
  
  n<-length(y)
  output<-matrix(ncol=3,nrow=size)
  for(i in 1:size)
  {
    alpha <- rnorm(1,mean(y)-beta*mean(x),1/(n*tau))
    m <- (sum(x*y)-alpha*n*mean(x))/sum(x**2)
    s <- 1/(tau*sum(x**2))
    beta <- rnorm(1,m,s)
    w <- y-alpha-beta*x
    tau <- rgamma(1,((n/2)+1),(sum(w**2)/2))
    output[i,] <- c(alpha,beta,tau)
  }
   data.frame(alpha=output[,1],beta=output[,2],tau=output[,3])
}


outLR<-LinReg(y=ob$y,x=ob$time,size=10000)
cat(paste0("alpha = ",format(round(mean(outLR$alpha),4),nsmall = 4), 
           " beta = ", format(round(mean(outLR$beta),4),nsmall = 4), 
           " tau = ",format(round(mean(outLR$tau),4),nsmall = 4)))
cat(paste0("sigma = ",format(round(sqrt(mean(1/outLR$tau)),4),nsmall = 4) ))

 จะเห็นได้ว่ามันดูยุ่งยากพอสมควรในการที่จะเขียน PDF หรือ likelihood สำหรับตัวแปรแต่ล่ะตัวครับ ถ้าเพียง 2-3 ตัวยังพอไหวแต่ถ้ามากกว่านี้ล่ะ คงเหนื่อยหน้าดูแถมในบางปัญหาก็ยากและไม่สามารถที่จะเขียน likelihood ได้ และนี่เองจึงเป็นที่มาของการพัฒนารูปแบบภาษาหรือคำสั่งในคอมพิวเตอร์ที่เรียกว่า Probabilistic Programming Language มาช่วยเรื่องพวกนี้ครับ

 

ต่อ