Metode Monte Carlo

Menghitung nili dari integral fungsi di atas dengan menggunakan metode Montecarlo dengan beberapa nilai m untuk dibandingkan mana nilai yang paling mendekati perhitungan dengan fungsi integrate dengan m = 10,100,1000,10000.

library(pracma)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

r1 <- function(x) {
  (3*exp(-2*x))}

curve(r1)

#Hitung nilai eksak dengan fungsi integrate
nilai.eksak<-integrate(r1,1,5)

#Metode Monte Carlo m = 10,100,1000,10000
mc_integral <- function(ftn, a, b,m=1000){
  #Membangkitkan x berdistribusi U(a,b)
  x <- runif(m,a,b)
  
  # Menghitung rata-rata dari output fungsi
  Gx <- ftn(x)
  Gx_m <- mean(Gx)
  theta.hat <- (b-a)*Gx_m
  return(theta.hat)
}

r1

## function(x) {
##   (3*exp(-2*x))}
## <bytecode: 0x0000000019b8a6f0>

set.seed(123)
int_r1 <- mc_integral(r1,a=1,b=5,m = 10000)
int_r2 <- mc_integral(r1,a=1,b=5,m = 1000)
int_r3 <- mc_integral(r1,a=1,b=5,m = 100)
int_r4 <- mc_integral(r1,a=1,b=5,m = 10)
int_r1

## [1] 0.2032603

int_r2

## [1] 0.1898719

int_r3

## [1] 0.1650022

int_r4

## [1] 0.2201995

nilai.eksak$value

## [1] 0.2029348

Untuk lebih memudahkan dalam membandingkan nilai-nilai tersebut, selanjutnya akan disajikan dalam tabel komparasi.

Komparasi Hasil

Berikut ini adalah hasil komparasi beberapa nilai m dengan nilai eksak yang dihitung dengan fungsi integrate.

hasil <-matrix (c(int_r1,int_r2,int_r3,int_r4, nilai.eksak$value), nrow =5, ncol =1)
rownames ( hasil ) <- c("m = 10000 ","m = 1000","m = 100","m = 10", "Nilai Eksak")
colnames ( hasil ) <- c("Nilai")
hasil<-as.data.frame(hasil)
hasil

Berdasarkan komparasi di atas, dapat disimpulkan bahwa semakin besar nilai mmaka akan semakin mendekati hasil hitungnya terhadap nilai eksak nya.

Metode Sweep Operator

Prosedur Sweep Operator

Secara numerik pendugaan koefisien regresi dilakukan dengan suatu algoritma Sweep Operator, Algoritma ini berdasarkan eliminasi (Gauss elimination) terhadap matriks simetrik A yang berukuran m×m. Unsur diagonal utama matriks A, yaitu \(a_{kk}≠0\) , k=1,2,…,m.

Prosedur ini digunakan untuk menentukan matriks kebalikan (inverse) \(A^{-1}\), penduga koefisien regresi \(\hat{B}\), dan JKG (Jumlah Kuadrat Galat).

Proses sweep dapat dilakukan secara bertahap untuk setiap model dan reversible.

Algoritma Sweep Operator

1. D = \(a_{kk}\).

2. Bagi baris ke-k dengan D.

3. Untuk setiap baris ke \(i ≠ k\), lakukan:

   1. Tetapkan \(B = a_{i,k}\).
   2. Kurangkan baris ke-i dengan B kali baris ke-k.
   3. \(a_{ik} = -B/D\).

4. \(a_{kk} = 1/D\).

Tahapan:

1. Menentukan matriks A

\[ A = \begin{bmatrix} 6 & 12 & 0 & 12\\ 12 & 28 & 0 & 23\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

Menggunakan algoritma sweep operator, maka

Tahapan pertama mencari \(y=b_0\) :

1. D = \(a_{kk}=a_{11}=6\).

\[ A = \begin{bmatrix} 6 & 12 & 0 & 12\\ 12 & 28 & 0 & 23\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

2. Bagi baris ke-1 dengan 6.

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ 12 & 28 & 0 & 23\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

3. Untuk setiap baris ke \(i ≠ k\), lakukan:

• Untuk \(i=2\)

  1. Tetapkan \(B = a_{ik}=a_{21}=12\).

  2. Kurangkan baris ke-2 dengan \(B=12\) kali baris ke-1.

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ 0 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{12}{6}=-2\).

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

• Untuk \(i=3\)

  1. Tetapkan \(B = a_{ik}=a_{31}=0\).

  2. Kurangkan baris ke-3 dengan \(B=0\) kali baris ke-1.

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{0}{6}=0\).

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ 12 & 23 & -2 & 28\ \end{bmatrix}\]

• Untuk \(i=4\)

  1. Tetapkan \(B = a_{ik}=a_{41}=12\).

  2. Kurangkan baris ke-4 dengan \(B=12\) kali baris ke-1.

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ 0 & -1 & -2 & 4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{12}{6}=-2\).

\[ A = \begin{bmatrix} 1 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

4. \(a_{kk} = \frac{1}{D}=\frac{1}{6}\).

\[ A = \begin{bmatrix} 1/6 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

Tahapan kedua mencari \(y=b_0+b_1X_1\) :

1. D = \(a_{kk}=a_{22}=4\).

\[ A = \begin{bmatrix} 1/6 & 2 & 0 & 2\\ -2 & 4 & 0 & -1\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

2. Bagi baris ke-2 dengan 4.

\[ A = \begin{bmatrix} 1/6 & 2 & 0 & 2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

3. Untuk setiap baris ke \(i ≠ k\), lakukan:

• Untuk \(i=1\)

  1. Tetapkan \(B = a_{ik}=a_{12}=2\).

  2. Kurangkan baris ke-1 dengan \(B=2\) kali baris ke-2.

\[ A = \begin{bmatrix} 7/6 & 0 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{2}{4}=-\frac{1}{2}\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

• Untuk \(i=3\)

  1. Tetapkan \(B = a_{ik}=a_{32}=0\).

  2. Kurangkan baris ke-3 dengan \(B=0\) kali baris ke-2.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{0}{4}=0\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -2 & -1 & -2 & 4\ \end{bmatrix}\]

• Untuk \(i=4\)

  1. Tetapkan \(B = a_{ik}=a_{42}=-1\).

  2. Kurangkan baris ke-4 dengan \(B=-1\) kali baris ke-21.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -5/2 & 0 & -2 & 15/4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=\frac{1}{4}\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -5/2 & 1/4 & -2 & 5/4\ \end{bmatrix}\]

4. \(a_{kk} = \frac{1}{D}=\frac{1}{4}\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

Tahapan ketiga mencari \(y=b_0+b_1X_1+b_2X_2\) :

1. D = \(a_{kk}=a_{33}=6\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 6 & -2\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

2. Bagi baris ke-3 dengan 6.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

3. Untuk setiap baris ke \(i ≠ k\), lakukan:

• Untuk \(i=1\)

  1. Tetapkan \(B = a_{ik}=a_{13}=0\).

  2. Kurangkan baris ke-1 dengan \(B=0\) kali baris ke-3.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{0}{6}=0\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

• Untuk \(i=2\)

  1. Tetapkan \(B = a_{ik}=a_{23}=0\).

  2. Kurangkan baris ke-2 dengan \(B=0\) kali baris ke-3.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=-\frac{0}{6}=0\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & -2 & 15/4\ \end{bmatrix}\]

• Untuk \(i=4\)

  1. Tetapkan \(B = a_{ik}=a_{43}=-2\).

  2. Kurangkan baris ke-4 dengan \(B=-2\) kali baris ke-3.

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & 0 & 37/12\ \end{bmatrix}\]

3. \(a_{ik} = -\frac{B}{D}=\frac{1}{3}\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1 & -1/3\\ -5/2 & 1/4 & 1/3 & 37/12\ \end{bmatrix}\]

4. \(a_{kk} = \frac{1}{D}=\frac{1}{6}\).

\[ A = \begin{bmatrix} 7/6 & -1/2 & 0 & 5/2\\ -1/2 & 1/4 & 0 & -1/4\\ 0 & 0 & 1/6 & -1/3\\ -5/2 & 1/4 & 1/3 & 37/12\ \end{bmatrix}\]

Selanjutnya akan dilakukan perhitungan koefisien regresi dari matriks berikut ini.

library(fastmatrix)

## 
## Attaching package: 'fastmatrix'

## The following objects are masked from 'package:pracma':
## 
##     geomean, hadamard, lu

library(microbenchmark)

y = c(1, 3, 3, 2, 2, 1) # vector y
x1 = c(3, 2, 1, 3, 2, 1)    # vector x1
x2 = c(-1, -1, -1, 1, 1, 1)  # vector x2
s = c(1,1,1,1,1,1)  # vector ‘satu’
X=cbind(s,x1,x2)    # menggabungkan vector s, x1, dan x2
XX=t(X)%*%X # X’X
Xy=t(X)%*%y # X’y
yX=t(y)%*%X # y’X
yy=t(y)%*%y # y’y
b1=cbind(XX,Xy) # menggabungkan X’X dan X’y
b2=cbind(yX,yy) # menggabungkan y’X dan y’y
b1=cbind(XX,Xy) # menggabungkan X’X dan X’y
b2=cbind(yX,yy) # menggabungkan y’X dan y’y
b=rbind(b1,b2)
b

##     s x1 x2   
## s   6 12  0 12
## x1 12 28  0 23
## x2  0  0  6 -2
##    12 23 -2 28

koef<-round(sweep.operator(b, 1:3),2)
coef_sweep <- koef[-4,4]
koef

##        s    x1   x2      
## s   1.17 -0.50 0.00  2.50
## x1 -0.50  0.25 0.00 -0.25
## x2  0.00  0.00 0.17 -0.33
##    -2.50  0.25 0.33  3.08

coef_sweep

##     s    x1    x2 
##  2.50 -0.25 -0.33

Maka didapatkan koefisien regresi b0, b1, dan b2 secara berturut-turut 2.5, -0.25, dan -0.33.

Tabel Koefisien Regresi

Selanjutnya akan dibuat tabel tahapan sweep operator. Untuk membuat tabel Tahapan SWEEP Operator, kita perlu mengatur nilai k.

Jika k=1 maka kita hanya mendapatkan dugaan intercept b0, k=2 maka akan mendapatkan b0 dan b1, jika k=3 maka akan mendapatkan b0,b1, dan b2.

library(fastmatrix)
library(microbenchmark)

coef_sweep_b0 <- sweep.operator(b,k = 1)
coef_sweep_b0

##             s x1 x2   
## s   0.1666667  2  0  2
## x1 -2.0000000  4  0 -1
## x2  0.0000000  0  6 -2
##    -2.0000000 -1 -2  4

#Koefisien b0
coef_sweep_b0[1,4]

## [1] 2

#SSE
coef_sweep_b0[4,4]

## [1] 4

coef_sweep_b1 <- sweep.operator(b,k = 1:2)
coef_sweep_b1

##            s    x1 x2      
## s   1.166667 -0.50  0  2.50
## x1 -0.500000  0.25  0 -0.25
## x2  0.000000  0.00  6 -2.00
##    -2.500000  0.25 -2  3.75

#Koefisien b0 & b1
coef_sweep_b1[1:2,4]

##     s    x1 
##  2.50 -0.25

#SSE
coef_sweep_b1[4,4]

## [1] 3.75

coef_sweep_b1[1,4]

## [1] 2.5

coef_sweep_b2 <- sweep.operator(b,k = 1:3)
coef_sweep_b2

##            s    x1        x2           
## s   1.166667 -0.50 0.0000000  2.5000000
## x1 -0.500000  0.25 0.0000000 -0.2500000
## x2  0.000000  0.00 0.1666667 -0.3333333
##    -2.500000  0.25 0.3333333  3.0833333

#Koefisien b0, b1, & b2
coef_sweep_b2[1:3,4]

##          s         x1         x2 
##  2.5000000 -0.2500000 -0.3333333

#SSE
coef_sweep_b2[4,4]

## [1] 3.083333

coef_sweep_b1_exc <- sweep.operator(b,k = c(1,3))
coef_sweep_b1_exc

##             s x1        x2           
## s   0.1666667  2 0.0000000  2.0000000
## x1 -2.0000000  4 0.0000000 -1.0000000
## x2  0.0000000  0 0.1666667 -0.3333333
##    -2.0000000 -1 0.3333333  3.3333333

#Koefisien b0 & b2
coef_sweep_b1_exc[c(1,3),4]

##          s         x2 
##  2.0000000 -0.3333333

#SSE
coef_sweep_b1_exc[4,4]

## [1] 3.333333

#Tabel Koefisien Regresi


coefb0_reg<-matrix(c(coef_sweep_b0[1,4],coef_sweep_b1[1,4],coef_sweep_b2[1,4],coef_sweep_b1_exc[1,4]), nrow = 4, ncol = 1)
colnames(coefb0_reg)<-c("Keofisien b0")
coefb0_reg

##      Keofisien b0
## [1,]          2.0
## [2,]          2.5
## [3,]          2.5
## [4,]          2.0

coefb1_reg<-matrix(c(0,coef_sweep_b1[2,4],coef_sweep_b2[2,4],0), nrow = 4, ncol = 1)
colnames(coefb1_reg)<-c("Keofisien b1")
coefb1_reg

##      Keofisien b1
## [1,]         0.00
## [2,]        -0.25
## [3,]        -0.25
## [4,]         0.00

coefb2_reg<-matrix(c(0,0,round(coef_sweep_b2[3,4],2),round(coef_sweep_b1_exc[3,4],2)), nrow = 4, ncol = 1)
colnames(coefb2_reg)<-c("Keofisien b2")
coefb2_reg

##      Keofisien b2
## [1,]         0.00
## [2,]         0.00
## [3,]        -0.33
## [4,]        -0.33

SSE_reg<-matrix(c(coef_sweep_b0[4,4],coef_sweep_b1[4,4],round(coef_sweep_b2[4,4],2),round(coef_sweep_b1_exc[4,4],2)), nrow = 4, ncol = 1)
colnames(SSE_reg)<-c("SSE")
SSE_reg

##       SSE
## [1,] 4.00
## [2,] 3.75
## [3,] 3.08
## [4,] 3.33

Step_reg<-matrix(c("Langkah1","Langkah2","Langkah3","Langkah4"), nrow = 4, ncol = 1)
colnames(Step_reg)<-c("Step")
Step_reg

##      Step      
## [1,] "Langkah1"
## [2,] "Langkah2"
## [3,] "Langkah3"
## [4,] "Langkah4"

tabel_koef<-cbind(Step_reg,coefb0_reg,coefb1_reg,coefb2_reg, SSE_reg)
tabel_koef<- as.data.frame(tabel_koef)
tabel_koef

Metode EM

Metode Em merupakan suatu prosedur iteratif untuk pendugaan parameter menggunakan metode kemungkinan maksimum (maximum likehood).

Konsep dasar: menggabungkan data pengamatan dengan data tidak teramati (hilang/tersembunyi) sehingga bentuk fungsi kemungkinannya lebih mudah untuk dianalisis.

Algoritma EM:

1. begin

2. initialize \(\theta^{0}, T , i = 0\)

3. do \(i = i+1\)

4. E step: compute \(Q(\theta;\theta^{i})\)

5. M step : \(\theta^{i+1}\) <- arg max \(Q(\theta;\theta^{i})\)

6. until \([Q(\theta^{i+1};\theta^{i})-Q(\theta^{i};\theta^{i-1})]\le T\)

7. return \(\hat{\theta} = \theta^{i+1}\)

8. end

Sebelum membuat syntax di R, terlebih dahlu perlu didefinisikan persamaan dari \(x_{1}d_{2}\), \(x_{2}d_{2}\) serta nilai maksimum dari fungsi \(Q(\theta)\) sebagai berikut:

E-Step

\[Q(\theta)=E(x_{2}|y,\theta^{0})log(\theta) + E(x_{3}|y,\theta^{0})log(1-\theta)\]

Pendugaan \(x_{1}\) dan \(x_{2}\) diperoleh sebagai berikut :

\[\hat{x}_{1}^{p} = 40*\frac{(\frac{1}{2})}{(\frac{1}{2}+\frac{\theta}{2})}\]

\[\hat{x}_{2}^{p} = 40*\frac{(\frac{\theta}{2})}{(\frac{1}{2}+\frac{\theta}{2})}\]

M-Step

Selanjutnya mendefinisikan fungsi yang memaksimumkan fungsi \(Q(\theta)\) sedemikian sehingga diperolah hasil sebagai berikut:

\[\theta^{1} = \frac{\hat{x}_{2}}{\hat{x}_{2}+x_{3}}\]

Setelah didapatkan persamaan-persamaan tersebut selanjutnya kita lanjutkan dengan syntax di R untuk menentukan nilai \(x_{1}\) dan \(x_{2}\) yang tersembunyi sebagai berikut:

#Definisikan nilai awal
x3<- 18
x2d1<- 1        
x2d2<- 1    
theta0<- 2
err<- 10
iterasi<- 1
while(err>10^-4){
  theta1<- (x2d1)/(x2d1+x3)
  x1d2<- 40*(0.5/(0.5+(theta1/2)))  
  x2d2<- 40*((theta1/2)/(0.5+(theta1/2)))   
  err<- abs(theta1-theta0)
  theta0<- theta1
  x2d1<-x2d2
  iterasi<- iterasi + 1
}

x1d2  

## [1] 29.00154

x2d2   

## [1] 10.99846

theta1

## [1] 0.3792373

iterasi

## [1] 15

Berdasarkan formula tersebut diperoleh nilai estimasi untuk \(x_{1}=29.0\) dan estimasi untuk nilai \(x_{2}=10.99\). Sementara itu, nilai estimasi \(\hat{\theta}=0.379\).