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 มาช่วยเรื่องพวกนี้ครับ

ต่อ

Probabilistic Programming Language

Like this:

Related

Share this post!

Like this:

Related