Integral Numerik dengan R

Fungsi Tiga Metode

Mendefinisikan fungsi untuk 3 metode yaitu Trapezoid, Simpson, dan Gauss Quadrature

1. Fungsi Trapezoid

trapezoid <- function(ftn, a, b, n, toleransi = 0.0001, ...){
  result <- data.frame(n=NULL, result=NULL, error=NULL)
  e <- 1 
  
  # proses memperoleh nilai exact
  exact <- integrate(ftn, a, b, ...)
  while(e > toleransi){
    
    # proses perhitungan trapezoid
    h <- (b-a)/n
    x.vec <- seq(a, b, by=h)
    f.vec <- sapply(x.vec, ftn, ...)
    Z <- h*(f.vec[1]/2 + sum(f.vec[2:n]) + f.vec[n+1]/2)
    
    # menampilkan nilai akhir
    e <- abs(Z-exact$value)
    result <- rbind(result, data.frame(n=n, result=Z, error=e))
    n <- n+1
  }
  result 
}

2. Fungsi Simpson

simpson <- function(ftn, a, b, n, toleransi = 0.0001, ...){
  result <- data.frame(n=NULL, result=NULL, error=NULL)
  e <- 1 
  
  # proses memperoleh nilai exact
  exact <- integrate(ftn, a, b, ...)
  while(e > toleransi){
    
    # proses perhitungan simpson
    n_s <- max(c(2*(n %/% 2), 4))
    h <- (b-a)/n_s
    
    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))
    
    # menampilkan nilai akhir
    e <- abs(S-exact$value)
    result <- rbind(result, data.frame(n=n, result=S, error=e))
    n <- n+1
  }
  
  result
}

3. Fungsi Gaussian Quadrature

gauss_quadrature <- function(ftn, a, b, n, toleransi = 0.0001, ...){
  result <- data.frame(n=NULL, result=NULL, error=NULL)
  e <- 1 
  
  # proses memperoleh nilai exact
  exact <- integrate(ftn, a, b, ...)
  while(e > toleransi){
    
    # proses mendapatkan koefisien dan titik gauss
    eval(parse(text=paste0("gL <- gaussLegendre(n, a=", a, ",", b, ")")))
    
    C_i <- gL$w # koefisien
    x_i <- gL$x # gauss point
    I <- sum(C_i * ftn(x_i, ...))
    
    # menampilkan nilai akhir
    e <- abs(I-exact$value)
    result <- rbind(result, data.frame(n=n, result=I, error=e))
    n <- n+1
  }
  
  result
}

4. Fungsi Monte Carlo Integration

monte_Carlo_integral <- function(ftn, a, b, m=1000, ...){
  # proses menghasilkan x berdistribusi Uniform(a,b)
  x <- runif(m, a, b)
  
  # menghitung rataan dari f(x)
  Gx <- ftn(x, ...)
  Gx_bar <- mean(Gx)
  
  theta.hat <- (b-a) * Gx_bar
  
  theta.hat
}

1. Distribusi Kumaraswamy

Diketahui pdf (probability density function) atau fkp (fungsi kepekatan peluang) dari distribusi Kumaraswamy adalah sebagai berikut \(f(x|a,b)=abx^{\alpha-1}(1-x^{(\alpha)})^{b-1}\), dimana \(x \in [0,1]\)

Bila diketahui \(a=5\) dan \(b=3\),

  1. Hitunglah \(E(X)\) dengan menggunakan metode Trapezoide (n=4), Simpson (n=4) dan four-point Gauss Quandrature.

  2. Bila diketahui nilai exaxct dari \(E(X)=\frac{b \Gamma (1+\frac{1}{a}) \Gamma(b)} {\Gamma (1+ \frac{1}{a}+b)}\), metode mana yang peling mendekati untuk menghitung nilai \(E(X)\) dengan toleransi 0.0001?

Penyelesaian

Subsitusi nilai a dan b ke dalam probability density function distribusi Kumaraswamy

\[f(x;5,3) = 5\ \cdot \ 3x^{5-1}(1-x^{5})^{3-1} \\ \ \ \ \ \ \ = 15x^{4} \cdot (1-x^{5})^{2}\]

a. Menghitung Nilai E(X) Secara Manual

a.1 Metode Trapezoid dengan n=4

Nilai harapan atau E(X) dari distribusi Kumaraswamy diperoleh dengan mengintegralkan pdf tersebut dengan interval \([0,1]\)

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

Menggunakan metode trapezoid dengan \(n=4\) diperoleh panjang setiap sub interval yaitu:

\[h = \frac{b-a}{n} \\ \ \ = \frac{1-0}{4} \\ \ \ \ \ \ \ \ \ =\frac{1}{4} = 0.25\] Sehingga bentuk dari metode trapezoidal adalah

\[T_{4} = \frac{h}{2} \ [E(x_0)\ +\ 2E(x_1)\ +\ 2E{x_2}\ +\ 2E(x_3)\ +\ 2E(x_4)]\] Selanjutnya, menghitung fungsi \(f(x)\) pada setiap \(x_i\)

\(E(x_0) = E(0) = 15(0)^5 \ (1-(0)^5)^2 = 0\)

\(E(x_1) = E(0.25) = 15(0.25)^5 \ (1-(0.25)^5)^2 = 0.01461984\)

\(E(x_2) = E(0.5) = 15(0.5)^5 \ (1-(0.5)^5)^2 = 0.43991089\)

\(E(x_3) = E(0.75) = 15(0.75)^5 \ (1-(0.75)^5)^2 = 2.07061679\)

\(E(x_4) = E(1) = 15(1)^5 \ (1-(1)^5)^2 = 0\)


Karena nilai \(h=0.25\) sehingga diperoleh \[T_{4} = \frac{h}{2} \ [E(x_0)\ +\ 2E(x_1)\ +\ 2E{x_2}\ +\ 2E(x_3)\ +\ E(x_4)] \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \frac{0.25}{2} \ [0 + 2(0.014621984) + 2(0.43991089) + 2(2.07061679) + 0] = 0.63128688\] sehingga

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

a.2 Metode Simpson dengan \(n=4\)

Oleh karena diketahui nilai \(n=4\), serta fungsi \(f(x_i)=\) dan nilai \(h\) diketahui, sehingga dari bentuk Simpson selanjutnya dihitung dan diperoleh nilai \(S_4\) sebagi berikut

\(\int_{0}^{1}x\ \cdot 15x^4(1-x^5)^2\ dx \approx S_4=\frac{h}{3}[E(x_0) + 4E(x_1) + 2E(x_2) + 4E(x_3) + 2E(x_4)]\)

\(S_4=\frac{h}{3}[E(x_0) + 4E(x_1) + 2E(x_2) + 4E(x_3) + E(x_4)] \\ \ \ \ \ = \frac{0.25}{3} [0 + 4(0.01461984) + 2(0.43991089) + 4(2.07061679) + 0] \\ \ \ \ \ = 0.76839736\)

sehingga

\(E(X) = \int_{0}^{1}x\ \cdot 15x^4(1-x^5)^2\ dx \\ \ \ \ \ \ \ \ \ \ = \int_{0}^{1} 15x^5(1-x^5)^2\ dx \approx\ S_4 = 0.76839736\)

a.3 Metode Four-Point Gaussian Quadrature

Akan dihitung nilai \(E(X)\) menggunakan metode four-point Gaussian Quadrature dengan bentuk umumnya adalah

\[\int_{a}^{b} f(x) dx \approx I = \sum_{i=1}^{n}C_if(x_i)\] dimana koefisien \(C_i\) adalah bobot dengan \(x_i\) adalah titik Gauss yang berbeda dengan interval \([-1,1]\), sehingga domain integral \([a,b]\) perlu ditransformasi dan mengubah variabel

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

\[dx = \frac{1}{2}(b-a)dt\] sehingga menjadi

\[ \int_{a}^{b} f(x)dx \to \int_{-1}^{1} f(t)dt = \int_{-1}^{1} f\left( \frac{1}{2} [t(b-a)+a+b] \right) \cdot \frac{1}{2}(b-a) dt\] Langkah 1

Untuk menghitung nilai harapan distribusi Kumaraswamy, mensubstitusikan nilai \(a\) dan \(b\) ke persamaan \(x\) dan \(dx\).

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

\(dx = \frac{1}{2}(1-0)dt = \frac{1}{2}dt\) sehingga diperoleh integral dari nilai harapan sebagai berikut:

\[E(X) = \int_{0}^{1} x \cdot 15x^4(1-x^5)^2\ dx \\ \ \ \ \ = \int_{0}^{1} 15x^5(1-x^5)^2\ dx \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \int_{-1}^{1} \frac{15}{2} \left( \frac{t+1}{2} \right)^5\ \left(1- \left( \frac{t+1}{2} \right)^5 \right)^2 \ dt \ \approx I\]

Langkah 2

Menghitung integral menggunakan Four-Point Gauss Quadrature dengan mengacu pada tabel gaussian quadrature untuk n=4

library(pracma)
gaussan <- gaussLegendre(n=4, a=-1, b=1)

daftar <- data.frame(i=c(1:4), gaussan[[2]], gaussan[[1]])
colnames(daftar) <- c("i", "Koefisien C_i", "Titik Gauss x_i")
daftar
##   i Koefisien C_i Titik Gauss 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

Selanjutnya mensubstitusikan nilai koefisien \(C_{i}\) dan titik gauss \(x_{i}\) ke dalam integral \(E(X)\) yang telah ditransformasi.

\(I = \int_{-1}^{1} E(t)dt = C_{1}E(t_{1}) + C_{2}E(t_{2}) + C_{3}E(t_{3}) + C_{4}E(t_{4}) \\ \ \ = C_{1}f(t_{1}) + C_{2}f(t_{2}) + C_{3}f(t_{3}) + C_{4}f(t_{4}) \\ \ \ = 0.3478548 \cdot f(-08611363) + 0.6521452 \cdot f(-0.3399810) + 0.6521452 \cdot f(0.3399810) + 0.3478548 \cdot f(0.8611363)\)

nilai fungsi \(f(x_i)\) diperoleh sebegai berikut

\(f(-0.8611363) = E(-0.8611363) = \frac{15}{2} \left( \frac{-0.8611363 + 1}{2} \right)^5\ \left(1- \left( \frac{-0.8611363 + 1}{2} \right)^5 \right)^2 \ = 0.0000121019\)

\(f(-0.3399810) = E(-0.3399810) = \frac{15}{2} \left( \frac{-0.3399810 + 1}{2} \right)^5\ \left(1- \left( \frac{-0.3399810 + 1}{2} \right)^5 \right)^2 \ = 0.029126407\)

\(f(0.3399810) = E(0.3399810) = \frac{15}{2} \left( \frac{0.3399810 + 1}{2} \right)^5\ \left(1- \left( \frac{0.3399810 + 1}{2} \right)^5 \right)^2 \ = 0.757589241\)

\(f(0.8611363) = E(0.8611363) = \frac{15}{2} \left( \frac{0.8611363 + 1}{2} \right)^5\ \left(1- \left( \frac{0.8611363 + 1}{2} \right)^5 \right)^2 \ = 0.477908858\)

sehingga nilai diperoleh nilai \(I\)

$I= 0.3478548 (0.0000121019) + 0.6521452 (0.029126407) + 0.6521452 (0.757589241) + 0.3478548 (0.477908858) \ = 0.679299934 $

Dengan semikian diperoleh

\(E(X) = \int_{-1}^{1} 15x^5(1-x^5)^2 \ dx \\ \ \ \ \ \ \ \ \ \ \ = \int_{-1}^{1} \frac{15}{2} \left( \frac{t+1}{2} \right)^5\ \left(1- \left( \frac{t+1}{2} \right)^5 \right)^2 \ dt \approx I = 0.679299934\)

a.4 Perbandingan hasil perhitungan manual dengan 3 metode

nilai_EX <- rbind(0.63128688, 0.76839736, 0.679299934)
nama_metode <- c('Trapezoid n=4', 'Simpson n=4', 'Four-point Gauss Quadrature')

data.frame(nama_metode, nilai_EX)
##                   nama_metode  nilai_EX
## 1               Trapezoid n=4 0.6312869
## 2                 Simpson n=4 0.7683974
## 3 Four-point Gauss Quadrature 0.6792999

b. Menghitung nilai E(X) dengan Nilai Exact

Nilai exact dari \(E(X)=\frac{b \Gamma (1+\frac{1}{a}) \Gamma(b)} {\Gamma (1+ \frac{1}{a}+b)}\).
Dari 2 metode (Trapezoid, Simpson, dan Gauss Quadrature) akan dicari metode mana yang hasilnya paling mendekati nilai exact \(E(X)\) dengan toleransi 0.0001?

Nilai exact distribusi Kumaraswamy adalah \(E(X) = \frac{b \Gamma (1+\frac{1}{a}) \Gamma(b)} {\Gamma (1+ \frac{1}{a}+b)} \\ \ \ \ \ \ \ \ \ \ \ = \frac{3 \cdot \Gamma (1+\frac{1}{5}) \Gamma(3)} {\Gamma (1+ \frac{1}{5}+3)} = 0.7102273\)

perhitungan (manual) menggunakan R

a <- 5
b <- 3
E_X <- (b*gamma(1+(1/a))*gamma(b)) / (gamma(1+(1/a)+b))
E_X
## [1] 0.7102273

Perbandingan nilai E(X) Nilai di atas dibandingkan dengan nilai yang dihasilkan dengan 3 metode.

E_Kumaraswamy <- function(x){15*(x^5) * (1-(x^5))^2}
batas_a <- 0  # batas bawah
batas_b <- 1  # batas atas
n <- 4

# Perhitungan menggunakan fungsi 3 metode yang telah disusun sebelumnya
kum_trapezoid <- trapezoid(E_Kumaraswamy, batas_a, batas_b, n)
kum_simpson <- simpson(E_Kumaraswamy, batas_a, batas_b, n)
kum_gauss <- gauss_quadrature(E_Kumaraswamy, batas_a, batas_b, n)

# Perbandingan hasil
compare_metode <- cbind(
  method = c('Trapezoid n=4', 'Simpson n=4', 'Four-Point Gauss Quadrature'),
  rbind(kum_trapezoid[which(kum_trapezoid$n==4),],
        kum_simpson[which(kum_simpson$n==4),],
        kum_gauss[which(kum_gauss$n==4),],
        make.row.names = FALSE)
)

# selisih hasil E(X) dengan perhitungan manual
compare_metode$selisih <- abs(compare_metode$result - (b*gamma(1+(1/a))*gamma(b)) / (gamma(1+(1/a)+b)))

# daftar perbandingan
compare_metode
##                        method n    result      error    selisih
## 1               Trapezoid n=4 4 0.6312869 0.07894039 0.07894039
## 2                 Simpson n=4 4 0.7683974 0.05817008 0.05817008
## 3 Four-Point Gauss Quadrature 4 0.6792999 0.03092735 0.03092735

Dari perbandingan nilai \(E(X\) di atas, dengan \(n=4\) dan toleransi = 0.0001 diperoleh bahwa selisih nilai \[E(X)\] yang dihasilkan oleh metode Guas Qaudrature paling kecil. Ini berarti bahwa metode tersebut merupakan metode yang paling mendekati nilai exact \(E(X)\) distribusi Kumaraswamy dengan nilai \(0.6792999\)

2. Distribusi 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]\]

Diketahui \(\lambda=2\) dan \(x=4\), hitunglah cdf(cumulative density function) dari distribusi Eksponensial di atas menggunakan metode Trapezoida, Simpson, Gauss Quadrature, dan Monte Carlo dan bandingkan mana yang terbaik.

Penyelesaian

CDF dari suatu distribusi didefinisikan sebagai berikut \[F(x) = \int_{0}^{\infty} f(x) dx\]

CDF dari distribusi Eksponensial \(x \sim exp(\lambda=2)\) adalah:

\[F(x) = \int_{0}^{\infty} \lambda e^{-\lambda t} \ dt\]

\[F(x=4|2) = \int_{0}^{4} 2e^{-2t} \ dt \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \left[ -e^{-2\ t} \ |_{0}^{4} \right] \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = \left[-e^{-8} - e^{0} \right] \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = 0,999664537\] Nilai exact juga dapat diperoleh dengan cara lain, yaitu sebagai berikut:

\[\ \ \ \ \ \ F(x) = 1-e^{- \lambda \ x} \\ F(x=4) = 1-e^{-2 \ \cdot \ 4} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ = 0,999664537\] Perhitungan nilai exact menggunakan R

# Cara 1
batas_a <- 0
batas_b <- 4
lambda <- 2
F_exp <- function(x, lambda) 1-exp(-lambda*x )
F_exp(batas_b, lambda)
## [1] 0.9996645
# Cara 2
f_exp <- function(x, lambda){
  lambda*exp(-1*lambda*x)
  }

# cdf fungsi distribusi Eksponensial
integrate(f_exp, batas_a, batas_b, lambda)
## 0.9996645 with absolute error < 1.1e-14

Perbandingan nilai E(X)

Selanjutnya dihitung CDF tersebut dengan menggunakan 4 metode yang disebutkan sebelumnya.

Pada simulasi ini, perbandingan metode integral numerik untuk distribusi Eksponensial, agar mendekati nilai exact, beberapa parameter yang digunakan:

  • Trapezoid dan Simpson menggunakan \(n=80\)

  • Gauss Quadrature menggunakan \(n=100\)

  • Monte Carlo dengan \(m = 100\)

  • \(Nilai \ toleransi = 0.00001\)

n_trap <- 80
n_simp <- 80
n_gauss <- 100
m_monte <- 100
toleransi <- 0.00001

F_exp_trapezoid <- trapezoid(f_exp, a=batas_a, b=batas_b, lambda, n=n_trap, toleransi=toleransi)

F_exp_simpson <- simpson(f_exp, a=batas_a, b=batas_b, lambda, n=n_simp, toleransi=toleransi)

F_exp_gauss <- gauss_quadrature(f_exp, a=batas_a, b=batas_b, lambda, n=n_gauss, toleransi=toleransi)

F_exp_monte <- monte_Carlo_integral(f_exp, a=batas_a, b=batas_b, lambda, m=m_monte)


# mengambil nilai E(X) tiap metode dengan parameter yg telah ditentukan sebelumnya
hasil <- data.frame(hasil= rbind(tail(F_exp_trapezoid, n=1)$result, 
                                         tail(F_exp_simpson, n=1)$result, 
                                         tail(F_exp_gauss, n=1)$result, 
                                         F_exp_monte))

# tabel perbandingan
perbandingan <- cbind(method=c('Trapezoid', 'Simpson', 'Gauss Quadrature', 'Monte Carlo'), hasil)

# menghilangkan kolom index
row.names(perbandingan) <- NULL

# menambahkan kolom 'selisih' -> selisih nilai E(X) yg dihasilkan tiap metode dengan perhitungan manual
perbandingan$selisih <- abs(perbandingan$hasil - F_exp(batas_b, lambda))

# menampilkan tabel akhir
perbandingan
##             method     hasil      selisih
## 1        Trapezoid 0.9996745 9.977400e-06
## 2          Simpson 0.9996651 5.547087e-07
## 3 Gauss Quadrature 0.9996645 6.994405e-15
## 4      Monte Carlo 0.9191832 8.048138e-02

Dari tabel di atas dapat dilihat bahwa selisih (error) terkecil diperoleh pada metode Gauss Quadrature. Dengan kata lain bahwa metode Gauss Quadrature n=100 dengan toleransi=0.00001 hasilnya paling mendekati nilai exact \(x \sim exp(\lambda = 2)\) atau \(F(x)= 0,999664537\)

Referensi