1 Acceptance Rejection (AR)

a. Acceptance Rejection Method Distribusi Poisson

PDF distribusi Poisson : \(P\left( x \right) = \frac{{e^{ - \lambda } \lambda ^x }}{{x!}}\) Nilai \(tx\)=\(x\)

Langkah-langkah untuk melakukan acceptance rejection adalah:

  1. Memilih nilai awal x=x0

  2. Generate \(\lambda\) interval antara \(0<\lambda<b\).

  3. generate \(u\) berdistribusi U(0,1)

  4. Jika \(u<f(x)/tx\) maka terima nilai \(\lambda\), else tolak \(\lambda\)

a=0
b=12
x=5
n=1500
lambda = runif(n,a,b)
U= runif(n,0,1)
fx=c()
f_t=c()
hasil=c()
fxx=c()
for (i in 1:n){
  fx[i]=((exp(-lambda[i]))*(lambda[i]^x)/factorial(x))
  fxx[i]=fx[i]/max(fx)
  if (U[i]<fxx[i]){
    hasil[i]="accept lambda"}
  else {hasil[i]="reject lambda"}
}
result = data.frame(lambda,U,fxx,hasil)
acc=subset(result,hasil== "accept lambda")

plot(lambda,U,pch=19)
points(lambda, fxx, col="red", lwd=3)
points(acc[,1], acc[,2], col="blue", lwd=3)
legend("topright", legend=c("reject","accept" ,"fx/tx"), col=c("black", "blue","red"), lty=1:1, cex=0.8,bg="white",lwd=3)

lambda_acc <- (as.numeric(acc[,1]))

Goodness of Fit untuk parameter Poisson

Menguji goodness of fit menggunakan package::fitdistrplus untuk mendapatkan estimasi parameter dan distribusi parameter.

#goodness of fit
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
descdist(lambda_acc, discrete = FALSE)

## summary statistics
## ------
## min:  1.329307   max:  11.76779 
## median:  5.69561 
## mean:  5.943938 
## estimated sd:  2.239328 
## estimated skewness:  0.3961287 
## estimated kurtosis:  2.604119
fit_w3  <- fitdist(lambda_acc, "weibull")
fit_g3  <- fitdist(lambda_acc, "gamma")
fit_nrm3 <- fitdist(lambda_acc,"norm")
test=gofstat(list(fit_w3, fit_g3, fit_nrm3), fitnames = c("weibull", "gamma", "normal")) 
test
## Goodness-of-fit statistics
##                                 weibull      gamma     normal
## Kolmogorov-Smirnov statistic 0.04071192 0.02996883 0.05442588
## Cramer-von Mises statistic   0.23563037 0.07669824 0.48565980
## Anderson-Darling statistic   1.60987871 0.69409096 3.27276893
## 
## Goodness-of-fit criteria
##                                 weibull   gamma   normal
## Akaike's Information Criterion 3137.067 3134.92 3162.662
## Bayesian Information Criterion 3146.198 3144.05 3171.792
test$aic
##  weibull    gamma   normal 
## 3137.067 3134.920 3162.662
test$bic
##  weibull    gamma   normal 
## 3146.198 3144.050 3171.792

Berdasarkan nilai AIC dan BIC minimum pada distribusi Gamma. maka \(\lambda\) ~ \(gamma(\alpha,\beta)\) dengan estimasi parameter berikut:

fit_g  <- fitdist(lambda_acc, "gamma")
summary(fit_g)
## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##       estimate Std. Error
## shape 6.607645 0.34220750
## rate  1.111664 0.05981655
## Loglikelihood:  -1565.46   AIC:  3134.92   BIC:  3144.05 
## Correlation matrix:
##           shape      rate
## shape 1.0000000 0.9624861
## rate  0.9624861 1.0000000
fit_g$estimate
##    shape     rate 
## 6.607645 1.111664

b. Acceptance Rejection Distribusi Lognormal

PDF Lognormal adalah \(f(x) = \frac{1}{(\sigma\sqrt{2 \pi})} e^{-(\frac{(log(x) - \mu)^2}{2 \sigma^2})}\).

Lognormal memiliki dua parameter yaitu \(f(x)\)~\(lognormal(\mu,\sigma)\)

Maka AR dilakukan pada masing-masing parameter.

AR parameter miu

Langkah-langkah untuk melakukan acceptance rejection adalah:

  1. Memilih nilai awal x=x0

  2. Generate \(\mu\) dengan interval antara \(a\) sampai dengan \(b\).

  3. generate \(u\) berdistribusi U(0,1)

  4. Jika \(u<f(x)/tx\) maka terima nilai \(\mu\), else tolak \(\mu\)

x=5
sigma0=3
a=-15
b=15
n=1000
mu0=runif(n,a,b)
u=runif(n,0,1)
tx=(0.1)*sqrt(2)/(sigma0*sqrt(pi))
fx0=c()
hasil0=c()

for (i in 1:n)
{
  fx0[i]=(exp(-1/(2*sigma0^2)*(log(x)-mu0[i])^2)/(sigma0*x*sqrt(2*pi)))/tx
  if (u[i]<=fx0[i])
  {hasil0[i]="accept mu"}
  else {hasil0[i]="reject mu"}
}
result0=cbind(mu0,u,fx0,hasil0)
acc0=subset(result0,hasil0== "accept mu")

plot(mu0,u,pch=20)
points(mu0, fx0, col="red", lwd=3)
points(acc0[,1],acc0[,2], col="blue", lwd=3)
legend("topright", legend=c("reject","accept" ,"fx/tx"), col=c("black", "blue","red"), lty=1:1, cex=0.8,bg="white",lwd=3)

miu= as.numeric(acc0[,1])
miu_acc = mean(as.numeric(acc0[,1]))

Goodness of Fit untuk parameter mu Lognormal

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

## summary statistics
## ------
## min:  -6.295455   max:  9.608832 
## median:  1.43101 
## mean:  1.485614 
## estimated sd:  2.946692 
## estimated skewness:  -0.0126816 
## estimated kurtosis:  2.831898

Hasil gambar menunjukkan bahwa observasi miu sangat dekat dengan distribusi normal. Parameter miu memiliki nilai dari \(-\infty<\mu<\infty\) , sehingga hanya didekati dist. normal.

fit_nrm <- fitdist(miu,"norm")
gofstat(fit_nrm,fitnames="normal")
## Goodness-of-fit statistics
##                                  normal
## Kolmogorov-Smirnov statistic 0.03135208
## Cramer-von Mises statistic   0.02192559
## Anderson-Darling statistic   0.20591000
## 
## Goodness-of-fit criteria
##                                  normal
## Akaike's Information Criterion 1222.813
## Bayesian Information Criterion 1229.808
summary(fit_nrm)
## Fitting of the distribution ' norm ' by maximum likelihood 
## Parameters : 
##      estimate Std. Error
## mean 1.485614  0.1882557
## sd   2.940648  0.1331168
## Loglikelihood:  -609.4067   AIC:  1222.813   BIC:  1229.808 
## Correlation matrix:
##      mean sd
## mean    1  0
## sd      0  1

AR parameter sigma

Langkah-langkah untuk melakukan acceptance rejection adalah:

  1. Memilih nilai awal x=x0

  2. Generate \(\sigma\) dengan interval antara \(a\) sampai dengan \(b\).

  3. generate \(u\) berdistribusi U(0,1)

  4. Jika \(u<f(x)/tx\) maka terima nilai \(\sigma\), else tolak \(\sigma\)

sigma= runif(n,0,b)
tx1=abs(exp(-1/2)/(sqrt(2*pi)*(log(x)-miu_acc)) )
fx1=c()
hasil1=c()

for (i in 1:n)
{
  fx1[i]=(exp(-1/(2*sigma[i]^2)*(log(x)-miu_acc)^2)/(x*sqrt(2*pi*sigma[i]^2)))/tx1
  if (u[i]<=fx1[i])
  {hasil1[i]="accept sigma"}
  else {hasil1[i]="reject sigma"}
}
result=cbind(sigma,u,fx1,hasil1)
acc=subset(result,hasil1== "accept sigma")
plot(sigma,u,pch=20)
points(sigma, fx1, col="red", lwd=3)
points(acc[,1], acc[,2], col="blue", lwd=3)
legend("topright", legend=c("reject","accept" ,"fx/tx"), col=c("black", "blue","red"), lty=1:1, cex=0.8,bg="white")

sigma_acc=(as.numeric(acc[,1]))^2
length(sigma_acc)
## [1] 15

Goodness of Fit untuk parameter sigma Lognormal

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

## summary statistics
## ------
## min:  0.009101063   max:  178.175 
## median:  5.536004 
## mean:  28.08668 
## estimated sd:  49.69205 
## estimated skewness:  2.365759 
## estimated kurtosis:  8.821015

Gambar menunjukkan observasi sigma berada di area beta. Namun distribusi beta memiliki batas x \(0<x<1\) maka selain itu sigma akan dicobakan distribusi Weibull dan gamma dan dibandingkan menurut Goodness of fit.

fit_w1  <- fitdist(sigma_acc, "weibull")
fit_g1 <- fitdist(sigma_acc, "gamma")
test11=gofstat(list(fit_w1, fit_g1), fitnames = c("weibull", "gamma")) 
test11
## Goodness-of-fit statistics
##                                weibull      gamma
## Kolmogorov-Smirnov statistic 0.1358459 0.19002292
## Cramer-von Mises statistic   0.0425115 0.05851064
## Anderson-Darling statistic   0.2613319 0.32252687
## 
## Goodness-of-fit criteria
##                                 weibull    gamma
## Akaike's Information Criterion 109.4648 109.7646
## Bayesian Information Criterion 110.8809 111.1807

Apabila dilihat dari kriteria AIC dan BIC parameter sigma dapat berdistribusi Gamma atau Weibull. AIC dan BIC denga nilai minimum maka memiliki kriteria lebih baik. Hal tersebut yang memiliki AIC dan BIC minimum adalah Distribusi Gamma. Maka sigma berdistribusi Weibull. \(\sigma^2\)~ \(gamma(\alpha,\beta)\). Estimasi parameter sbb:

summary(fit_g1)
## Fitting of the distribution ' gamma ' by maximum likelihood 
## Parameters : 
##         estimate  Std. Error
## shape 0.30702479 0.088486560
## rate  0.01092747 0.005927279
## Loglikelihood:  -52.88231   AIC:  109.7646   BIC:  111.1807 
## Correlation matrix:
##           shape      rate
## shape 1.0000000 0.5238502
## rate  0.5238502 1.0000000

2. Slice Sampling

Slice sampling atau sampling terpotong adalah metode Markov Chain Monte Carlo yang digunakan untuk sampel point dari distribusi yang diberikan. Algoritma nya adalah :

  1. Menentukan interval awal \(a\) dan \(b\) pada sebuah distribusi.

  2. Melakukan random distribusi uniform secara horisontal antara interval \(a\) sampai \(b\)

  3. Didapatkan nilai random \(x\) maka didapatkan \(f(x)\)

  4. Mengambil irisan secara vertikal dengan interval \([x,f(x)]\)

  5. Generate titik sampel baru sehingga dapat nilai \(a*\) dan \(b*\) baru.

  6. Dilakukan seterusnya iterasi \(i+1\)

a. Slice Sampling Distribusi Laplace

pdf distribusi Laplace : \(f(x|\mu,B)\)= \(\frac{1}{2b}e^\frac{{-|x-\mu|}}{b}\).

Syarat \(|x|,x\) bila \(x>0\) dan \(-x\) bila \(x<0\).

n=1000 
mu=3 #parameter inisial awal
B=5 #parameter inisial awal

a = 2 #a dan b awal
b = 15
rand_hor=c(); rand_ver=c(); aa=c(); bb=c()

for (i in 1:n){
  rand_1 = runif(1, a, b)
  rand_hor[i]=rand_1
  fx = (exp(-abs(rand_1-mu)/B))/(2*B)
  rand_2 = runif(1,0,fx)
  rand_ver[i]=rand_2
  
  a=mu+B*log(rand_2*2*B) #untuk x<0
  b=mu-B*log(rand_2*2*B) #untuk x>0
  
  aa[i]= a
  bb[i]= b
}
hasil= cbind(rand_hor,rand_ver,aa,bb)
head(hasil)
##        rand_hor    rand_ver         aa        bb
## [1,]  8.6025625 0.016539900  -5.996973 11.996973
## [2,] -1.6098850 0.019426017  -5.192785 11.192785
## [3,]  5.5016903 0.005479962 -11.520360 17.520360
## [4,] -2.6140609 0.009851392  -8.587787 14.587787
## [5,]  7.5572114 0.010396902  -8.318311 14.318311
## [6,] -0.0523316 0.053303894  -0.145804  6.145804
plot(rand_hor,rand_ver,pch=20)

Selanjutnya mendapatkan distribusi dari random horisontal.

Menguji distribusi \(x\) berdasarkan uji Anderson Darling atau Kolmogorov Smirnov menggunakan package::lawstat. H0 : Data berdistribusi Laplace

H1 : Data tidak berdistribusi Laplace

library(lawstat)
## Warning: package 'lawstat' was built under R version 3.6.3
#uji Anderson Darling
p_val_AD=laplace.test(rand_hor)$A2
p_val_KS=laplace.test(rand_hor)$D
p_value=cbind(p_val_AD,p_val_KS)
p_value
##       p_val_AD  p_val_KS
## [1,] 0.7146043 0.7937733

b. Slice Sampling Distribusi Eksponensial

PDF exponensial adalah

\(f(x)=\lambda x^{-\lambda x}\),dimana \(x>0\)

n=1000 
lambda=2
#interval
a=0
b=10 #karena eksponensial nilai x>0
rand_hor=c(); rand_ver=c()
aa=c()
bb=c()

for (i in 1:n){
  rand_1= runif(1,a,b)
  rand_hor[i]=rand_1
  fx=lambda*exp(-lambda*rand_1)
  rand_2=runif(1,0,fx)
  rand_ver[i]=rand_2
  a=0
  b=-log(rand_2/lambda)/lambda
  aa[i]=a
  bb[i]=b
  
}
hasil=cbind(rand_hor,rand_ver,aa,bb)
plot(rand_hor,rand_ver,pch=20)

Menguji distribusi \(x\) berdasarkan uji Anderson Darling,Kolmogorov Smirnov dan Cramer–von Mises statistic dan kriteria perbandingan distribusi menggunakan AIC dan BIC.

#goodness of fit
library(fitdistrplus)
#win.graph()
descdist(rand_hor,discrete=FALSE)

## summary statistics
## ------
## min:  0.0005181801   max:  9.227536 
## median:  0.3413687 
## mean:  0.5034231 
## estimated sd:  0.5783308 
## estimated skewness:  4.989053 
## estimated kurtosis:  59.36753
fit_w  <- fitdist(rand_hor, "weibull")
fit_g  <- fitdist(rand_hor, "gamma")
fit_lg <- fitdist(rand_hor,"lnorm")
fit_exp <- fitdist(rand_hor,"exp")
test1=gofstat(list(fit_w, fit_g,fit_lg,fit_exp), fitnames = c("weibull", "gamma","lognormal","exponensial")) 
test1
## Goodness-of-fit statistics
##                                 weibull      gamma lognormal exponensial
## Kolmogorov-Smirnov statistic 0.01439527 0.01231342 0.0737675  0.01599120
## Cramer-von Mises statistic   0.02571092 0.03909230 1.6242438  0.06328975
## Anderson-Darling statistic   0.32562750 0.32202747 9.8992272  0.40473723
## 
## Goodness-of-fit criteria
##                                 weibull    gamma lognormal exponensial
## Akaike's Information Criterion 629.9193 631.1068  753.0685    629.3514
## Bayesian Information Criterion 639.7348 640.9223  762.8840    634.2592
#hasil dari slice sampling hasil random horisontal sebagai nilai x bisa berdistribusi weibull atau exponential. 
#Karena setiap running berbeda keputusan.