1. Slice Sampling Regresi Sederhana satu Data

Diketahui sebuah persamaan regresi \(y=\beta_{0}+\beta_{1}x+\epsilon\) dimana \(\epsilon\)~\(N(0,\sigma^2)\). Dilakukan simulasi Monte Carlo Markov Chain (MCMC) dengan metode slice sampling untuk bisa mendapatkan konvergen distribusi dari parameter \(\beta_{0}\), \(\beta_{1}\). pdf untuk \(\epsilon\)

\(f(\epsilon) = \frac{1}{(\sigma\sqrt{2 \pi})} e^{-(\frac{(y - (\beta_{0}+\beta_{1}x))^2}{2 \sigma^2})}\) Slice sampling ini mengasumsikan bahwa \(\sigma^2\) konstan.

n=100#k1
y=15 #k2 
x=3 #k3

#inisial value
b0 = 1 #k30
b1= 1 #40
sigma =1 #k50


q=-log(0.6283185307*10^(-3)*sigma^2)
a_b0=((-b1*x)+y-(sigma*sqrt(q)))
b_b0=((-b1*x)+y+(sigma*sqrt(q))) 
a_b1=(-b0+y-(sigma*sqrt(q)))/x 
b_b1=(-b0+y+(sigma*sqrt(q)))/x 


randhor_b1=c(); randver_b1=c()
randhor_b0=c(); randver_b0=c()
ba_b1=c(); bb_b1=c(); ba_b0=c(); bb_b0=c()

for (i in 1:n){
  random11 = runif(1,a_b0,b_b0)
  b0=random11
  randhor_b0[i]=random11
  fxx =exp(-0.5*((y-b0-(b1*x))/sigma)^2)/(sqrt(2*pi*sigma^2))
  random22= runif(1,0,fxx)
  randver_b0[i]=random22
  
  q=-log(0.6283185307*10^(-3)*sigma^2)
  a_b0=((-b1*x)+y-(sigma*sqrt(q)))
  b_b0=((-b1*x)+y+(sigma*sqrt(q))) 
  
  ba_b0[i]=b_b0
  bb_b0[i]=a_b0
  
  random1= runif(1,a_b1,b_b1)
  b1=random1
  randhor_b1[i]=random1
  fx=exp(-0.5*((y-b0-b1*x)/sigma)^2)/(sqrt(2*pi*sigma^2))
  random2 = runif(1,0,fx)
  randver_b1[i]=random2
  
  q=-log(0.6283185307*10^(-3)*sigma^2)
  a_b1=(-b0+y-(sigma*sqrt(q)))/x 
  b_b1=(-b0+y+(sigma*sqrt(q)))/x 
  
  
  ba_b1[i]= b_b1
  bb_b1[i]= a_b1
  
}
result_b1= cbind(randhor_b1,randver_b1,bb_b1,ba_b1)
head(result_b1)
##      randhor_b1   randver_b1      bb_b1    ba_b1
## [1,]  4.6351007 1.233550e-21  0.5226220 2.332774
## [2,]  2.1920592 5.237333e-06 -0.1469671 1.663185
## [3,]  0.8261420 2.496266e-42  4.5072470 6.317399
## [4,]  5.2940533 1.689982e-29  0.5861522 2.396304
## [5,]  1.8634313 7.540402e-05 -0.3281621 1.481990
## [6,]  0.7613518 4.558621e-57  5.2030730 7.013225
result_b0=cbind(randhor_b0,randver_b0,bb_b0,ba_b0)
head(result_b0)
##      randhor_b0   randver_b0     bb_b0     ba_b0
## [1,]  10.716906 6.209577e-02  9.284772 14.715228
## [2,]  12.725673 4.106189e-31 -1.620530  3.809926
## [3,]  -1.236969 9.707702e-22  5.708594 11.139050
## [4,]  10.526315 7.595002e-03  9.806346 15.236802
## [5,]  13.269258 5.047934e-46 -3.597388  1.833068
## [6,]  -3.324447 8.256050e-37  6.694478 12.124934
#win.graph()
par(mfrow=(c(2,2)))
hist(randhor_b1)
hist(randhor_b0)
plot(randhor_b0,randhor_b1,pch=20)

library(fitdistrplus)
## Warning: package 'fitdistrplus' was built under R version 3.6.3
## Loading required package: MASS
## Loading required package: survival
## Loading required package: npsurv
## Loading required package: lsei
#win.graph()
descdist(randhor_b0,discrete=FALSE)

## summary statistics
## ------
## min:  -12.38487   max:  33.30951 
## median:  4.905349 
## mean:  8.013319 
## estimated sd:  12.12308 
## estimated skewness:  0.5731379 
## estimated kurtosis:  2.126702
#win.graph()
descdist(randhor_b1,discrete=FALSE)

## summary statistics
## ------
## min:  -6.69324   max:  8.486957 
## median:  3.615741 
## mean:  2.286277 
## estimated sd:  4.073341 
## estimated skewness:  -0.6271582 
## estimated kurtosis:  2.17506

Hasil pengujian Goodness of fit menggunakan package::fitdistplus, baik random horisontal \(\beta_{0}\) dan \(\beta_{1}\) menunjukkan bahwa data berdistribusi normal. Karena simbol observasi berada pada bintang.

cat("Uji normalitas untuk beta 1")
## Uji normalitas untuk beta 1
ks.test(randhor_b1,"pnorm",mean(randhor_b1),sd=sd(randhor_b1))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  randhor_b1
## D = 0.13442, p-value = 0.05391
## alternative hypothesis: two-sided
cat("Uji normalitas untuk beta 0")
## Uji normalitas untuk beta 0
ks.test(randhor_b0,"pnorm",mean(randhor_b0),sd=sd(randhor_b0))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  randhor_b0
## D = 0.12633, p-value = 0.08218
## alternative hypothesis: two-sided

p-value kedua parameter menunjukkan lebih dari 0,05 maka slice sampling untuk random horisontal \(\beta_{0}\) dan \(\beta_{1}\) konvergen dalam distribusi normal.

#plot untuk b0
fit_nrm <- fitdist(result_b0[,1],"norm")
#win.graph()
par(mfrow=(c(2,2)))
hist(randhor_b0)
denscomp(fit_nrm, legendtext = "normal") #density histogram
qqcomp  (fit_nrm, legendtext = "normal")
ppcomp  (fit_nrm, legendtext = "normal")

#plot untuk b1
fit_nrm_b1 <- fitdist(result_b1[,1],"norm")
#win.graph()
par(mfrow=c(2,2))
hist(randhor_b1)
denscomp(fit_nrm_b1, legendtext = "normal") #density histogram
qqcomp  (fit_nrm_b1, legendtext = "normal")
ppcomp  (fit_nrm_b1, legendtext = "normal")

cat("parameter konvergen distribusi dari b0")
## parameter konvergen distribusi dari b0
fit_nrm
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##       estimate Std. Error
## mean  8.013319  1.2062310
## sd   12.062310  0.8529341
cat("parameter konvergen distribusi dari b1")
## parameter konvergen distribusi dari b1
fit_nrm_b1
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters:
##      estimate Std. Error
## mean 2.286277  0.4052923
## sd   4.052923  0.2865849

Kesimpulan : apabila dicobakan n=1000 tidak bisa mendapatkan parameter yang konvergen distribusi normal. Sehingga hanya dicobakan n=100 atau paling besar n=150.

2. Slice Sampling untuk Regresi N data

Selanjutnya mendapatkan konvergen distribusi parameter dari N data regresi. Maka menggunakan fungsi likelihood sebagai berikut

\(L(\beta_{0},\beta_{1},\sigma^2)\)=\(\prod_{i=1}^{n}f(\epsilon_{i})\) \(L(\beta_{0},\beta_{1},\sigma^2)\)=\(\prod_{i=1}^{n}\frac{1}{(\sigma\sqrt{2 \pi})} e^{-(\frac{(y -(\beta_{0}+\beta_{1}x))^2}{2 \sigma^2})}\) \(L(\beta_{0},\beta_{1},\sigma^2)\)= \((\frac{1}{(\sigma\sqrt{2 \pi})})^n e^{-\sum_{i=1}^{n}(\frac{(y -(\beta_{0}+\beta_{1}x))^2}{2 \sigma^2})}\) Slice sampling ini mengasumsikan bahwa \(\sigma^2\) konstan.

#masukkan data
xn=c(2,3,4,2,5,3,4,5,3,4)
yn=c(5,6,7,8,5,6,6,5,9,4)
N=400
#inisial value parameter
b0n = 1 #k1
b1n= 1 #k2
sigman =1 #k3
nn=length(yn)

sum_y = sum(yn) #k6
sum_x =sum(xn) #k5
sum_y2 = sum(yn^2) #k8
sum_x2 = sum(xn^2) #k7
sum_xy =sum(xn*yn) #k9


a_b1n =-(1/sum_x2)*(((b0n*sum_x)-sum_xy+sqrt(abs((b0n^2*sum_x^2)-(2*b0n*sum_x*sum_xy)+sum_xy^2+(9.210340372*sum_x2*sigman^2)-(sum_x2*sigman^2*nn*log(6.283185307*sigman^2))-(sum_x2*nn*b0n^2)+(2*sum_x2*b0n*sum_y)-(sum_x2*sum_y2)))))
b_b1n =-(1/sum_x2)*(((b0n*sum_x)-sum_xy-sqrt(abs((b0n^2*sum_x^2)-(2*b0n*sum_x*sum_xy)+sum_xy^2+(9.210340372*sum_x2*sigman^2)-(sum_x2*sigman^2*nn*log(6.283185307*sigman^2))-(sum_x2*nn*b0n^2)+(2*sum_x2*b0n*sum_y)-(sum_x2*sum_y2)))))

a_b0n=-(1/nn)*((sum_x*b1n-sum_y+sqrt(abs((sum_x^2*b1n^2)-(2*b1n*sum_x*sum_y)+sum_y^2-(nn*b1n^2*sum_x2)+(9.210340372*nn*sigman^2)-(nn*sum_y2)+(2*nn*b1n*sum_xy)-(sigman^2*nn^2*log(6.283185307*sigman^2))))))
b_b0n=-(1/nn)*((sum_x*b1n-sum_y-sqrt(abs((sum_x^2*b1n^2)-(2*b1n*sum_x*sum_y)+sum_y^2-(nn*b1n^2*sum_x2)+(9.210340372*nn*sigman^2)-(nn*sum_y2)+(2*nn*b1n*sum_xy)-(sigman^2*nn^2*log(6.283185307*sigman^2))))))

randhor_b1n=c(); randver_b1n=c()
randhor_b0n=c(); randver_b0n=c()
ba_b1n=c(); bb_b1n=c(); ba_b0n=c(); bb_b0n=c()

for (i in 1:N){
  random11n = runif(1,a_b0n,b_b0n)
  b0n=random11n
  randhor_b0n[i]=random11n
  fxxn =(1/(sqrt(2*pi*sigman^2)))^nn*exp((-0.5*(nn*b0n^2/sigman^2))+sum((yn*b1n*xn/sigman^2)-(0.5*b1n^2*xn^2/sigman^2)-(b0n*b1n*xn/sigman^2)-(0.5*yn^2/sigman^2)+(yn*b0n/sigman^2)))
  random22n= runif(1,0,fxxn)
  randver_b0n[i]=random22n
  
a_b0n=-(1/nn)*((sum_x*b1n-sum_y+sqrt(abs((sum_x^2*b1n^2)-(2*b1n*sum_x*sum_y)+sum_y^2-(nn*b1n^2*sum_x2)+(9.210340372*nn*sigman^2)-(nn*sum_y2)+(2*nn*b1n*sum_xy)-(sigman^2*nn^2*log(6.283185307*sigman^2))))))
b_b0n=-(1/nn)*((sum_x*b1n-sum_y-sqrt(abs((sum_x^2*b1n^2)-(2*b1n*sum_x*sum_y)+sum_y^2-(nn*b1n^2*sum_x2)+(9.210340372*nn*sigman^2)-(nn*sum_y2)+(2*nn*b1n*sum_xy)-(sigman^2*nn^2*log(6.283185307*sigman^2))))))

  ba_b0n[i]=b_b0n
  bb_b0n[i]=a_b0n
  
  random1n= runif(1,a_b1n,b_b1n)
  b1n=random1n
  randhor_b1n[i]=random1n
  fxn=(1/(sqrt(2*pi*sigman^2)))^nn*exp((-0.5*(nn*b0n^2/sigman^2))+sum((yn*b1n*xn/sigman^2)-(0.5*b1n^2*xn^2/sigman^2)-(b0n*b1n*xn/sigman^2)-(0.5*yn^2/sigman^2)+(yn*b0n/sigman^2)))
  random2n = runif(1,0,fxn)
  randver_b1n[i]=random2n
  
 a_b1n =-(1/sum_x2)*(((b0n*sum_x)-sum_xy+sqrt(abs((b0n^2*sum_x^2)-(2*b0n*sum_x*sum_xy)+sum_xy^2+(9.210340372*sum_x2*sigman^2)-(sum_x2*sigman^2*nn*log(6.283185307*sigman^2))-(sum_x2*nn*b0n^2)+(2*sum_x2*b0n*sum_y)-(sum_x2*sum_y2)))))
b_b1n =-(1/sum_x2)*(((b0n*sum_x)-sum_xy-sqrt(abs((b0n^2*sum_x^2)-(2*b0n*sum_x*sum_xy)+sum_xy^2+(9.210340372*sum_x2*sigman^2)-(sum_x2*sigman^2*nn*log(6.283185307*sigman^2))-(sum_x2*nn*b0n^2)+(2*sum_x2*b0n*sum_y)-(sum_x2*sum_y2)))))
  
  
  ba_b1n[i]= b_b1n
  bb_b1n[i]= a_b1n
  
}
result_b1n= cbind(randhor_b1n,randver_b1n,bb_b1n,ba_b1n)
cat("hasil random horisontal dan vertikal b1 \n")
## hasil random horisontal dan vertikal b1
head(result_b1n)
##      randhor_b1n   randver_b1n       bb_b1n    ba_b1n
## [1,]   0.8252604  3.017604e-20  0.476191979 1.8428994
## [2,]   1.3142597  3.287963e-28  0.002218211 1.1261972
## [3,]   0.1361426  3.049874e-15 -0.055088418 1.0442238
## [4,]   0.5248633  6.382338e-44  0.699971362 2.2004136
## [5,]   2.1305451 7.679708e-100 -0.154217207 0.9052308
## [6,]   0.2420501  8.687275e-12 -0.411661177 0.5623021
result_b0n=cbind(randhor_b0n,randver_b0n,bb_b0n,ba_b0n)
cat("hasil random horisontal dan vertikal b0 \n")
## hasil random horisontal dan vertikal b0
head(result_b0n)
##      randhor_b0n   randver_b0n     bb_b0n   ba_b0n
## [1,]   1.5080121  1.896581e-17  0.2855145 4.914485
## [2,]   3.7702964  2.622610e-14  1.0218878 5.401289
## [3,]   4.0349284  3.431187e-31 -1.0551042 4.055286
## [4,]   0.4035542  1.491713e-69  3.8337274 7.413275
## [5,]   4.4873599  6.166685e-12  2.2683803 6.257576
## [6,]   5.6280680 3.421602e-132 -4.6038659 1.890050
par(mfrow=(c(2,2)))
hist(randhor_b1n)
hist(randhor_b0n)
plot(randhor_b0n,randhor_b1n,pch=20,col="blue")

library(fitdistrplus)
descdist(randhor_b0n,discrete=FALSE)

## summary statistics
## ------
## min:  -3.810668   max:  16.8816 
## median:  6.728157 
## mean:  6.648429 
## estimated sd:  3.448928 
## estimated skewness:  -0.1818341 
## estimated kurtosis:  2.891255
descdist(randhor_b1n,discrete=FALSE)

## summary statistics
## ------
## min:  -3.203052   max:  2.2867 
## median:  -0.1876152 
## mean:  -0.1684294 
## estimated sd:  0.9472947 
## estimated skewness:  0.1371414 
## estimated kurtosis:  3.044984
cat("Uji normalitas untuk beta 1")
## Uji normalitas untuk beta 1
ks.test(randhor_b1n,"pnorm",mean(randhor_b1n),sd=sd(randhor_b1n))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  randhor_b1n
## D = 0.026573, p-value = 0.9402
## alternative hypothesis: two-sided
cat("Uji normalitas untuk beta 0")
## Uji normalitas untuk beta 0
ks.test(randhor_b0n,"pnorm",mean(randhor_b0n),sd=sd(randhor_b0n))
## 
##  One-sample Kolmogorov-Smirnov test
## 
## data:  randhor_b0n
## D = 0.041597, p-value = 0.4932
## alternative hypothesis: two-sided
#plot untuk b0
fit_nrmn <- fitdist(result_b0n[,1],"norm")
#win.graph()
par(mfrow=(c(2,2)))
hist(randhor_b0n)
denscomp(fit_nrmn, legendtext = "normal") #density histogram
qqcomp  (fit_nrmn, legendtext = "normal")
ppcomp  (fit_nrmn, legendtext = "normal")

#plot untuk b1
fit_nrm_b1n <- fitdist(result_b1n[,1],"norm")
#win.graph()
par(mfrow=c(2,2))
hist(randhor_b1n)
denscomp(fit_nrm_b1n, legendtext = "normal") #density histogram
qqcomp  (fit_nrm_b1n, legendtext = "normal")
ppcomp  (fit_nrm_b1n, legendtext = "normal")