Metode Trapezoidal

Misal, \(T\) adalah nilai hampiran untuk \(E(x) = \int_a^b g(x)dx\) dengan metode Trapezoidal, maka

\(\begin{aligned} h &= \frac{b-a}{n}, \quad \mathrm{sehingga\ } a = x_0 < x_1< x_2 < ... < x_{n-1} < x_n = b \\ T &= \frac{h}{2}\sum_{i=0}^n cg(x_i), \quad \mathrm{di\ mana\ } c = 1 \mathrm{\ untuk\ } i \in (0,n) \mathrm{,\ selainnya\ } c = 2.\\ \end{aligned}\)

Nilai hampiran \(E(x)\) menggunakan metode Trapezoidal dengan \(n=4\) dapat dilakukan secara manual dengan tahapan berikut:

\(\begin{aligned} h &= \frac{b-a}{n} \\ &= \frac{1-0}{4} = 0.25, \quad \mathrm{sehingga\ } x_0 = 0,\ a_1 = 0.25,\ a_2 = 0.5,\ a_3 = 0.75,\ a_4 = 1.\\ \end{aligned}\)

Perhitungan berikutnya diringkas dalam tabel.

Sehingga, nilai hampiran dari \(E(x)\) adalah

\(\begin{aligned} T &= \frac{0.25}{2} (0 + 0.0292 + 0.8798 + 4.1412) \\ &= 0.125(5.0502) \\ &= 0.6313 \\ \end{aligned}\)

Metode Simpson

Misal, \(S\) adalah nilai hampiran untuk \(E(x) = \int_a^b g(x)dx\) dengan metode Simpson, maka

\(\begin{aligned} h &= \frac{b-a}{n}, \quad \mathrm{sehingga\ } a = x_0 < x_1< x_2 < ... < x_{n-1} < x_n = b \\ S &= \frac{h}{3}\sum_{i=0}^n cg(x_i), \quad \mathrm{di\ mana\ } c = 1 \mathrm{\ untuk\ } i \in (0,n)\mathrm{,\ } c = 4 \mathrm{\ untuk\ } i \mathrm{\ ganjil,\ } c = 2 \mathrm{\ untuk\ } i \mathrm{\ genap}. \\ \end{aligned}\)

Nilai hampiran \(E(x)\) menggunakan metode Simpson dengan \(n=4\) dapat dilakukan secara manual dengan tahapan berikut:

\(\begin{aligned} h &= \frac{b-a}{n} \\ &= \frac{1-0}{4} = 0.25, \quad \mathrm{sehingga\ } x_0 = 0,\ a_1 = 0.25,\ a_2 = 0.5,\ a_3 = 0.75,\ a_4 = 1.\\ \end{aligned}\)

Perhitungan berikutnya diringkas dalam tabel.

Sehingga, nilai hampiran untuk \(E(x)\) adalah

\(\begin{aligned} S &= \frac{0.25}{3} (0 + 0.0585 + 0.8798 + 8.2825) \\ &= 0.0833(9.2208) \\ &= 0.7684 \\ \end{aligned}\)

Metode Gauss (Adaptive) Quadrature

Misal, \(Q\) adalah nilai hampiran untuk \(E(x) = \int_a^b g(x)dx\) dengan metode Gauss (Adaptive) Quadrature, maka

\[ Q = \frac{b-a}{2}\sum_{i=1}^n w_{i}f\left(\frac{b-a}{2}x_{i}+\frac{b+a}{2}\right) \]

Untuk \(n = 4\), nilai \(w\) dan \(x\) adalah (Sumber: pomax.github.io)

Nilai hampiran \(E(x)\) menggunakan metode four-point Gauss (Adaptive) Quadrature dapat dilakukan secara manual seperti diringkas dalam tabel berikut:

di mana \(k_i = 0.5x_{i}+0.5\)

Sehingga, nilai hampiran bagi \(E(x)\) adalah

\(\begin{aligned} Q &= \frac{1-0}{2} (0.0380 + 0.9881 + 0 + 0.3325) \\ &= 0.5(1.3586) \\ &= 0.6793 \\ \end{aligned}\)

Nilai Harapan Exact

Nilai exact bagi nilai harapan distribusi Kumaraswamy, yaitu

\[ E(X) = \frac{b\Gamma(1+\frac{1}{a})\Gamma(b)}{\Gamma(1+\frac{1}{a}+b)} \]

dapat dihitung dengan menggunakan R dengan sintaks berikut:

ExKumaraswamy <- function(a,b){
  E <- (b * gamma(1 + (1/a)) * gamma(b)) / gamma(1 + (1/a) + b)
  return(E)
}

Untuk \(a=5\) dan \(b=3\), diperoleh nilai exact:

Ex <- ExKumaraswamy(a = 5, b = 3)
Ex

## [1] 0.7102273

Perbandingan nilai hampiran menggunakan teknik integrasi numerik dengan nilai exact ditampilkan dalam tabel berikut.

Berdasarkan tabel di atas, four-point Gauss (Adaptive) Quadrature memberikan nilai hampiran yang paling mendekati dibandingkan dua metode lainnya.

Fungsi R untuk Integrasi Numerik

Untuk menghitung nilai integral dengan teknik numerik dari suatu fungsi \(f\) pada interval \(a\) dan \(b\) dapat dilakukan dengan custom function berikut ini:

IntTrapezoid <- function(f, a, b, n){
  
  #' Metode Trapezoid, dengan banyaknya selang n
  
  h <- (b - a)/n
  c <- c(1, rep(2, n - 1), 1)
  x <- seq(a, b, by=h)
  
  Int <- (h/2) * sum(c * sapply(x, f))
  return(Int)
}


IntSimpson <- function(f, a, b, n){
  
  #' Metode Simpson, dengan banyaknya selang n
  
  h <- (b - a)/n
  c <- c(1, rep(c(4,2), length.out = n-1), 1)
  x <- seq(a, b, by=h)
  
  Int <- (h/3) * sum(c * sapply(x, f))
  return(Int)
}


library(gaussquad)
IntGaussQuad <- function(f, a, b, n){
  
  #' Metode Gauss (Adaptive) Quadrature, n points
  
  Lq <- legendre.quadrature.rules(max(2,n))[[n]]
  k <- ((b - a)/2)*Lq$x + ((b + a)/2)
  Int <- ((b - a)/2) * sum(Lq$w * f(k))
  return(Int)
}

Selanjutnya nilai harapan dari distribusi Kumaraswamy

\[ E(x) = \int_0^1 15x^5(1-x^5)^2 dx \] dapat dihitung dengan mengintegralkan fungsi

gx <- function(x) (15*x^5)*(1-x^5)^2

pada interval \(a=0\) dan \(b=1\).

Hasil integral dan selisih terhadap nilai exact adalah sebagai berikut:

ExNum <- c(
    `Trapezoidal, n = 4`        = IntTrapezoid(f = gx, a = 0, b = 1, n = 4),
    `Simpson, n = 4`            = IntSimpson(f = gx, a = 0, b = 1, n = 4),
    `Gauss Quadrature, 4 point` = IntGaussQuad(f = gx, a = 0, b = 1, n = 4) )

data.frame(Nilai.integral = ExNum,
           Nilai.exact    = Ex,
           Selisih        = abs(ExNum - Ex) )

Perhitungan R memberikan output yang relatif sama dengan hasil perhitungan manual.

Integrasi Numerik Iteratif

Perhitungan integral di atas menggunakan nilai \(n\) yang ditentukan. Selanjutnya, nilai integral dapat dihitung hingga toleransi tertentu (misalnya 0.0001) dengan cara mencari \(n\) optimal. Hal ini dapat dilakukan secara iteratif menggunakan fungsi berikut:

IterIntegral <- function(fun.method, f, a, b, 
                         exact, tol = 1e-4, max.iter = 1e3){
  
  #' Integrasi integrasi numerik
  #'
  #' fun.method : fungsi metode integrasi numerik
  #' f : fungsi yang akan diintegralkan pada interval a dan b
  #' exact : nilai exact
  #' tol : nilai toleransi
  #' max.iter: total iterasi maksimum
  
  err <- n <- 1

  while( err > tol & n <= max.iter){
    int <- fun.method(f, a, b, n)
    err <- abs(int - exact)
    n <- n+1
  }
  
  out <- c(int = int, exact = exact, err = err, n = n)
  return(out)
  
}

Hasil integrasi numerik secara iteratif adalah sebagai berikut:

result <- data.frame(rbind(
  IterIntegral(IntTrapezoid, f = gx, a = 0, b = 1, exact = Ex),
  IterIntegral(IntSimpson,   f = gx, a = 0, b = 1, exact = Ex),
  IterIntegral(IntGaussQuad, f = gx, a = 0, b = 1, exact = Ex) ))

names(result) <- c("Nilai.integral", "Nilai.exact", "Selisih", "Iterasi")
cbind(Metode = c("Trapezoidal", "Simpson", "Gauss Quadrature"), result)

Berdasarkan tabel di atas, ketiga metode menghasil nilai hampiran yang lebih kecil dari toleransi. Namun, metode Simpson dengan \(n = 35\) memberikan nilai hampiran yang paling mendekati nilai exact (dengan selisih 8.23e-5).

Nilai CDF Exact

Distribusi Eksponensial mempunyai CDF

\[ F(x) = 1 - \mathrm{\ exp}(-\lambda x) \] Dengan R, CDF dapat ditulis sebagai

CDFExp <- function(lambda, x){ 1 - exp(-lambda * x) }

Sehingga, nilai exact dari CDF sebuah distribusi Eksponensial dengan \(\lambda = 2\) dan \(x = 4\) adalah

Ex <- CDFExp(lambda = 2, x = 4)
Ex

## [1] 0.9996645

CDF dengan Integrasi Numerik

Fungsi distribusi kumulatif (CDF) dapat dihitung dengan cara mengintegralkan fungsi kepekatan peluang (PDF). Distribusi Eksponensial dengan \(\lambda = 2\) mempunyai PDF

\[ f(x\ |\ \lambda = 2) = 2 \mathrm{\ exp}(-2 x) \] Ditulis dalam bentuk fungsi R

PDFExp2 <- function(x){ 2 * exp(-2 * x)}

Sehingga untuk \(x = 4\), CDF adalah hasil integral pada interval \(0 \le x \le 4\).

Fungsi IntTrapezoid(), IntSimpson() dan IntGaussQuad() akan digunakan kembali untuk menghitung nilai integral teknik numerik dengan metode Trapezoidal, Simpson dan Gauss (Adaptive) Quadrature. Demikian juga dengan fungsi IterIntegral() untuk melakukan integrasi secara iteratif.

Untuk metode Monte Carlo, dibuat fungsi baru sebagai berikut

IntMonteCarlo <- function(f, a, b, n){
  x <- runif(n, a, b)
  gx <- f(x)
  gx_mean <- mean(gx)
  theta_hat <- (b-a) * gx_mean
  return(theta_hat)
}

Hasil integrasi numerik secara iteratif menggunakan metode Trapezoidal, Simpson, Gauss (Adaptive) Quadrature dan Monte Carlo serta selisih terhadap nilai exact yang diperoleh sebelumnya ditampilkan dalam tabel berikut:

set.seed(999)
result <- data.frame(rbind(
  IterIntegral(IntTrapezoid, f = PDFExp2, a = 0, b = 4, exact = Ex),
  IterIntegral(IntSimpson,   f = PDFExp2, a = 0, b = 4, exact = Ex),
  IterIntegral(IntGaussQuad, f = PDFExp2, a = 0, b = 4, exact = Ex),
  IterIntegral(IntMonteCarlo, f = PDFExp2, a = 0, b = 4, exact = Ex)))

names(result) <- c("Nilai.integral", "Nilai.exact", "Selisih", "Iterasi")
cbind(Metode = c("Trapezoidal", "Simpson", "Gauss Quadrature", "Monte Carlo"), 
      result)

Berdasarkan tabel di atas, keempat metode menghasil nilai hampiran yang lebih kecil dari toleransi, dan metode Monte Carlo dengan \(m = 963\) memberikan nilai hampiran yang paling mendekati nilai exact (dengan selisih 5.1e-6). Meskipun demikian, karena metode Monte Carlo bergantung pada pengacakan pada saat membangkitkan sebaran seragam, maka hasilnya bisa berubah tergantung dari random seed yang digunakan.

Hal menarik yang dapat dilihat dari tabel di atas adalah perbandingan antara metode Trapezoidal dan Simpson. Untuk kasus ini, metode Simpson memberikan hampiran nilai yang lebih baik (6.44e-5 : 9.99e-5) dan juga jumlah selang yang lebih sedikit dibandingkan Trapezoidal (22 : 232).