Tugas Praktikum Sains Data 10

1. Persamaan Kumaraswamy

Diketahui PDF (Probability Density Function) dari Kumaraswamy adalah sebagai berikut :

\(f(x|a,b) = abx^{a-1}(1-x^{a})^{b-1}\), dimana \(x\in(0,1)\)

Jika diketahui \(a = 5\) dan \(b = 3\)

a. Hitunglah \(E(X)\) dengan menggunakan metode Trapezoida untuk \(n = 4\), Simpson untuk \(n = 4\), dan four-point Gauss Quadrature tanpa menggunakan software apapun. (Boleh diketik ataupun ditulis tangan kemudian difoto)

b. Jika diketahui nilai exact dari \(E(X) = \frac{b\Gamma(1+\frac{1}{a})\Gamma(b)}{\Gamma(1+\frac{1}{a}+b)}\), metode mana yang paling mendekati untuk menghitung nilai \(E(X)\) dengan toleransi \(0.0001\)? Gunakan R!

Jawab :

a. Menghitung dengan manual

Diketahui \(a = 5\) dan \(b = 3\) maka,
\(f(x) = 15x^{4}(1-x^{5})^{2}\)

Maka untuk mencari \(E(X)\) dari distribusi Kumaraswamy adalah

\[E(X) = \int_{0}^{1} x \cdot 15x^{4}(1-x^{5})^{2} \; dx = \int_{0}^{1} 15x^{5}(1-x^{5})^{2} \; dx\]

Metode Trapezoida (\(n = 4\))

Diketahui fungsi untuk menghitung nilai ekspektasi dari distribusi Kumaraswamy :
\(g(x) = 15x^{5}(1-x^{5})^{2}\)
Panjang setiap sub interval adalah sebagai berikut :
\(h = \frac{b - a}{n} = \frac{1 - 0}{4} = \frac{1}{4} = 0.25\)
Karena \(n = 4\) maka bentuk dari metode trapezoidal adalah :
\(T_{4} = \frac{h}{2}[g(x_{0}) + 2 \cdot g(x_{1}) + 2 \cdot g(x_{2}) + 2 \cdot g(x_{3}) + g(x_{4})]\)
Menghitung fungsi \(g(x)\) pada setiap \(x_{i}\)
\(g(x_{0}) = g(0) = 0\)
\(g(x_{1}) = g(0.25) = 0.01461984\)
\(g(x_{2}) = g(0.5) = 0.4399109\)
\(g(x_{3}) = g(0.75) = 2.070617\)
\(g(x_{4}) = g(1) = 0\)
Karena \(h = 0.25\) diperoleh

\(T_{4} = \frac{0.25}{2}[0 + 2 \cdot 0.01461984 + 2 \cdot 0.4399109 + 2 \cdot 2.070617 + 0] = 0.6312869\)
Sehingga dengan menggunakan metode trapezoidal adalah :
\(E(X) = \int_{0}^{1} 15x^{5}(1-x^{5})^{2} \; dx = 0.6312869\)
Metode Simpson (\(n = 4\))

Diketahui fungsi untuk menghitung nilai ekspektasi dari distribusi Kumaraswamy :
\(g(x) = 15x^{5}(1-x^{5})^{2}\)
Panjang setiap sub interval adalah sebagai berikut :
\(h = \frac{b - a}{n} = \frac{1 - 0}{4} = \frac{1}{4} = 0.25\)
Karena \(n = 4\) maka bentuk dari metode simpson adalah :
\(S_{4} = \frac{h}{3}[g(x_{0}) + 4 \cdot g(x_{1}) + 2 \cdot g(x_{2}) + 4 \cdot g(x_{3}) + g(x_{4})]\)
Menghitung fungsi \(g(x)\) pada setiap \(x_{i}\)
\(g(x_{0}) = g(0) = 0\)
\(g(x_{1}) = g(0.25) = 0.01461984\)
\(g(x_{2}) = g(0.5) = 0.4399109\)
\(g(x_{3}) = g(0.75) = 2.070617\)
\(g(x_{4}) = g(1) = 0\)
Karena \(h = 0.25\) diperoleh

\(S_{4} = \frac{0.25}{3}[0 + 4 \cdot 0.01461984 + 2 \cdot 0.4399109 + 4 \cdot 2.070617 + 0] = 0.7683974\)
Sehingga dengan menggunakan metode simpson adalah :
\(E(X) = \int_{0}^{1} 15x^{5}(1-x^{5})^{2} \; dx = 0.7683974\)
Metode Four-Point Gauss Quadrature

Diketahui fungsi untuk menghitung nilai ekspektasi dari distribusi Kumaraswamy :
\(g(x) = 15x^{5}(1-x^{5})^{2}\)
Karena pada soal yang diinginkan adalah metode four-point gaussian quadratur, berarti nilai \(n = 4\). Berdasarkan tabel gaussian quadratur nilai-nilai koefisien dan titik gauss adalah sebagai berikut:

library(pracma)
fpgauss <- gaussLegendre(n=4,a=-1,1)
data.frame(i = c(1:4), C_i = fpgauss[[2]], x_i = fpgauss[[1]])

##   i       C_i        x_i
## 1 1 0.3478548 -0.8611363
## 2 2 0.6521452 -0.3399810
## 3 3 0.6521452  0.3399810
## 4 4 0.3478548  0.8611363

Sebelum kita menggunakan nilai-nilai koefisien dan titik gauss, perlu dilakukan transformasi sehingga domain integral menjadi \([−1,1]\).

\(x = \frac{1}{2}[t(b - a) + a + b]\) \(x = \frac{1}{2}[t(1 - 0) + 1 + 0]\) \(x = \frac{1}{2}[t + 1]\) dan

\(dx = \frac{1}{2}(b - a)dt\) \(dx = \frac{1}{2}(1 - 0)dt = \frac{1}{2}dt\)
Sehingga integralnya menjadi

\(\int_{0}^{1} 15x^{5}(1-x^{5})^{2} \; dx = \int_{-1}^{1} \frac{15}{2}(\frac{t+1}{2})^{5}[1 - (\frac{t+1}{2})^{5}]^{2} \; dt\)

Dari transformasi tersebut kita menemukan \(h(x) = \frac{15}{2}(\frac{t+1}{2})^{5}[1 - (\frac{t+1}{2})^{5}]^{2}\)

Kemudian, kita menggunakan nilai - nilai koefisien dan titik gauss saat \(n = 4\)

\(\int_{-1}^{1} \frac{15}{2}(\frac{t+1}{2})^{5}[1 - (\frac{t+1}{2})^{5}]^{2} \; dt \approx I = C_{1} \cdot h(t_{1}) + C_{2} \cdot h(t_{2}) + C_{3} \cdot h(t_{3}) + C_{4} \cdot h(t_{4})\)

\(I = 0.3478548 \cdot h(−0.8611363) + 0.6521452 \cdot h(−0.3399810) + 0.6521452 \cdot h(0.3399810) + 0.3478548 \cdot h(0.8611363)\)

Menghitung fungsi \(h(x)\) pada setiap \(x_{i}\)

\(h(−0.8611363) = \frac{15}{2}(\frac{−0.8611363+1}{2})^{5}[1 - (\frac{−0.8611363+1}{2})^{5}]^{2} = 1.21019 \cdot 10^{-5}\)
\(h(−0.3399810) = \frac{15}{2}(\frac{−0.3399810+1}{2})^{5}[1 - (\frac{−0.3399810+1}{2})^{5}]^{2} = 0.02912642\)
\(h(0.3399810) = \frac{15}{2}(\frac{0.3399810+1}{2})^{5}[1 - (\frac{0.3399810+1}{2})^{5}]^{2} = 0.7575892\)
\(h(0.8611363) = \frac{15}{2}(\frac{0.8611363+1}{2})^{5}[1 - (\frac{0.8611363+1}{2})^{5}]^{2} = 0.4779089\)
sehingga,

\(I = 0.3478548 \cdot 1.21019 \cdot 10^{-5} + 0.6521452 \cdot 0.02912642 + 0.6521452 \cdot 0.7575892 + 0.3478548 \cdot 0.4779089 = 0.6792999\)
Sehingga dengan menggunakan metode four-point gauss quadrature adalah :
\(E(X) = \int_{0}^{1} 15x^{5}(1-x^{5})^{2} \; dx = 0.6792999\)

b. Mengetahui metode terbaik untuk menghitung E(X)

Metode Trapezoida

trapezoid <- function(ftn, a, b, tol) {
  hasil_T <- data.frame(n=NULL, hasil=NULL, error=NULL)
  err = 1
  n = 4
  eksak <- integrate(ftn, a, b)
  while (err>tol){
  h <- (b-a)/n
  x.vec <- seq(a, b, by = h)
  f.vec <- sapply(x.vec, ftn)     # ftn(x.vec)
  Trap <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
  err <- abs(Trap - eksak$value)
  hasil_T <- rbind(hasil_T, data.frame(n=n, hasil=Trap, error=err))
  n = n + 1
  }
  hasil_T
}

Metode Simpson

simpson_n <- function(ftn, a, b, tol) {
  hasil_S <- data.frame(n=NULL, hasil=NULL, error=NULL)
  err = 1
  n = 4
  eksak <- integrate(ftn, a, b)
  while (err>tol){
  nb <- max(c(2*(n %/% 2), 4))
  h <- (b-a)/nb
  x.vec1 <- seq(a+h, b-h, by = 2*h)       # ganjil
  x.vec2 <- seq(a+2*h, b-2*h, by = 2*h)   # genap
  f.vec1 <- sapply(x.vec1, ftn)   # ganjil
  f.vec2 <- sapply(x.vec2, ftn)   # genap
  S <- h/3*(ftn(a) + ftn(b) + 4*sum(f.vec1) + 2*sum(f.vec2))
  err <- abs(S - eksak$value)
  hasil_S <- rbind(hasil_S, data.frame(n=n, hasil=S, error=err))
  n = n+1
  }
  hasil_S
}

Metode Gauss Quadrature

gaussian_q <- function(ftn, a, b, tol) {
  hasil_G <- data.frame(n=NULL, hasil=NULL, error=NULL)
  err = 1
  n = 2
  eksak <- integrate(ftn, a, b)
  while (err>tol){
    Lq <- gaussLegendre(n, a, b)
    xi = Lq$x
    wi = Lq$w 
    G <- sum(wi * ftn(xi))
    err <- abs(G - eksak$value)
    hasil_G <- rbind(hasil_G, data.frame(n=n, hasil=G, error=err))
    n = n + 1
  }
  hasil_G
}

Perbandingan hasil ketiga metode

f <- function(x){
  15*x^5*(1-x^5)^2
}
a <- 0
b <- 1
Tol <- 0.0001

Trapez <- trapezoid(f, a, b, Tol)
Simps <- simpson_n(f, a, b, Tol)
Gauss <- gaussian_q(f, a, b, Tol)

Perbandingan <- cbind(
  Metode = c('Trapezoidal', 'Simpson', 'Gauss Quadrature'), 
  rbind(tail(Trapez, n=1), tail(Simps, n=1), tail(Gauss, n=1), make.row.names = F)
)
Perbandingan

##             Metode  n     hasil        error
## 1      Trapezoidal 23 0.7101278 9.950361e-05
## 2          Simpson 34 0.7103096 8.227772e-05
## 3 Gauss Quadrature  6 0.7101422 8.504999e-05

Sehingga jika kita lihat dari hasil ketiga metode, metode Simpson mempunyai error terkecil. Sehingga metode terbaik untuk menghitung ekspektasi dari distribusi Kumaraswamy adalah metode Simpson

2. Persamaan Eksponensial

Diketahui PDF (Probability Density Function) dari distribusi eksponensial adalah sebagai berikut :

\(f(x|a,b) = \lambda e^{-\lambda x}\), dimana \(x\in(0,\infty)\)

Jika diketahui \(\lambda = 2\). Hitunglah CDF dari distribusi eksponensial tersebut untuk \(x = 4\), dengan menggunakan metode Trapezoida, Simpson, Gauss Quadrature dan Monte Carlo dengan R. Metode mana yang paling baik? Silahkan tentukan nilai toleransi, \(n\) dan \(m\) sendiri.

Jawab :

Diketahui \(\lambda = 2\), maka \(f(x) = 2e^{-2x}, x\in(0,\infty)\)

CDF dari distribusi eksponensial adalah :
\(F(x) = \int_{0}^{x} 2e^{-2x} \; dx\)
Akan dicari CDF distribusi eksponensial untuk \(x = 4\)
\(F(x) = \int_{0}^{4} 2e^{-2x} \; dx = [-e^{-2x}]^{4}_{0} = (-e^{-2 \cdot 4}) - (-e^{-2 \cdot 0}) = 0.9996645\)

atau dengan R :

fe <- function(x){
  2*exp(-2*x)
}
a <- 0
b <- 4
eksak_E <- integrate(fe, a, b)
eksak_E

## 0.9996645 with absolute error < 1.1e-14

Akan dibandingkan hasil \(F(x) = \int_{0}^{4} 2e^{-2x} \; dx\) menggunakan 4 metode, yaitu Trapezoidal, Simpson, Gauss Quadrature, dan Monte Carlo. Untuk metode Trapezoidal, Simpson, dan Gauss Quadrature dapat menggunakan fungsi yang telah dibuat pada nomor 1.

Metode Monte Carlo

mc_integral <- function(ftn, a, b, m){
  #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)
}

Perbandingan hasil keempat metode

Untuk perbandingan metode integral untuk distribusi eksponensial akan digunakan nilai toleransi = \(10^{-6}\) dan contoh acak sebanyak \(m = 10^{6}\). Untuk nilai \(n\) akan mengikuti nilai toleransi, sehingga akan didapatkan nilai \(n\) yang maksimum agar hasilnya mendekati nilai eksak dari distribusi eksponensial.

Tol <- 1e-6
m <- 1e6
Trapez_E <- trapezoid(fe, a, b, Tol)
Simps_E <- simpson_n(fe, a, b, Tol)
Gauss_E <- gaussian_q(fe, a, b, Tol)
Monte_E <- mc_integral(fe, a, b, m)
Perbandingan <- cbind(
  Metode = c('Trapezoidal', 'Simpson', 'Gauss Quadrature', 'Monte Carlo'),
  rbind(tail(Trapez_E, n=1), tail(Simps_E, n=1), tail(Gauss_E, n=1), c("-",Monte_E, abs(Monte_E-eksak_E$value)),make.row.names = F)
)
Perbandingan

##             Metode    n             hasil                error
## 1      Trapezoidal 2310 0.999665536518128 9.99146030578046e-07
## 2          Simpson   70 0.999665483336256 9.45964158338697e-07
## 3 Gauss Quadrature    7  0.99966448380689 5.35652073807569e-08
## 4      Monte Carlo    -  1.00257805764377  0.00291352027167424

Sehingga jika kita lihat dari hasil keempat metode, metode Gauss Quadrature mempunyai error terkecil. Sehingga metode terbaik untuk menghitung CDF dari distribusi eksponensial adalah metode Gauss Quadrature