Numerical Integration
Soal
Diketahui pdf (probability density function) dari distirbusi Kumaraswamy adalah sebagai berikut \[f(x|a,b)=ab\ x^{a-1}(1-x^a)^{b-1}\ ,\ dimana \ x \in (0,1)\] Jika diketahui a = 5 dan b = 3
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).
Jika diketahui nilai exact dari \(E(x)= \frac{b\ \Gamma(1+\frac{1}{2})\ \Gamma(b)}{\Gamma(1+\frac{1}{a}+b)}\), metode mana yang paling mendekati untuk menghitung nilai \(E(x)\) dengan toleransi 0.0001? Gunakan R!
Diketahui pdf (probability density function) dari distirbusi Eksponensial adalah sebagai berikut \[f(x|\lambda)=\lambda \ exp(-\lambda x)\ ,\ dimana \ x \in (0,\infty)\] Jika diketahui \(\lambda =2\). Hitunglah CDF dari distribusi exponensial 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.
Jawaban Soal 1
Pdf dari distirbusi Kumaraswamy adalah \[f(x|a,b)=ab\ x^{a-1}(1-x^a)^{b-1}\] Karena \(a=5\) dan \(b =3\), maka fungsinya dapat ditulis:
\[f(x)=15x^4(1-x^5)^2\] dengan batas \(0<x<1\). Sehingga untuk \(E(x)\) dapat diperoleh sebagai berikut: \[E(x)=\int_{0}^{1}x\cdot f(x)\ dx=\int_{0}^{1} 15x^5(1-x^5)^2\ dx\] Dengan kata lain, akan dicari luas daerah di bawah kurva \(f(x)=15x^5(1-x^5)^2\)
x <- seq(0,1,0.01)
y <- x*(15*x^4 * (1-x^5)^2)
plot(x,y,type = "l", ylab = "y",main=expression(f(x)==15*x^5 * (1-x^5)^2))
abline(h=0,v=c(0,1))
polygon(x,y,col = "grey")
a. Menghitung E(x) secara manual
Metode Trapezoida
Jadi, melalui metode Trapezoida untuk \(n=4\), diperoleh \(E(x)\approx0.6312869\)
Metode Simpson
Jadi, melalui metode Simpson untuk \(n=4\), diperoleh \(E(x)\approx0.7683974\)
Metode Four-Point Gauss Quadrature
Jadi, melalui metode Four-Point Gauss Quadrature, diperoleh \(E(x)\approx0.6792999\)
Tabel Hasil Perhitungan Manual
| Metode | n | Hasil |
|---|---|---|
| Trapezoida | 4 | 0.6312869 |
| Simpson | 4 | 0.7683974 |
| Gauss Quadrature | 4 | 0.6792999 |
b. Menghitung E(x) Menggunakan R
Pada bagian ini akan dihitung \(E(x)\) menggunakan metode Trapezoida, Simpson, dan Gauss Quadrature dengan toleransi 0.0001.
Diketahui bahwa Mean dari distirbusi Kumaraswamy adalah \(E(x)= \frac{b\ \Gamma(1+\frac{1}{2})\ \Gamma(b)}{\Gamma(1+\frac{1}{a}+b)}\)
a=5;b=3
meank = b * gamma(1+1/a) * gamma(b) / gamma(1+1/a+b)
meank## [1] 0.7102273Jadi, secara exact, \(E(x)=0.7102273\)
Metode Trapezoida
trapezoid <- function(ftn, a, b, n = 100) {
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)
return(Trap)
}f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
while(err>tol){
res_trap <- trapezoid(f,0,1,n = n)
err <- abs(res_trap-exact_value)
cat("n=",n,", result=",res_trap,", error=",err,"\n",sep = "")
n=n+1
if(n==1000){break}
}## n=4, result=0.6312869, error=0.07894042
## n=5, result=0.6738123, error=0.03641497
## n=6, result=0.6914928, error=0.01873454
## n=7, result=0.6997126, error=0.01051473
## n=8, result=0.7039057, error=0.006321603
## n=9, result=0.7062116, error=0.004015656
## n=10, result=0.7075597, error=0.002667593
## n=11, result=0.7083885, error=0.001838822
## n=12, result=0.7089199, error=0.001307439
## n=13, result=0.7092729, error=0.0009544148
## n=14, result=0.7095147, error=0.000712649
## n=15, result=0.7096846, error=0.0005426739
## n=16, result=0.7098069, error=0.0004204042
## n=17, result=0.7098966, error=0.0003306616
## n=18, result=0.7099637, error=0.0002636072
## n=19, result=0.7100146, error=0.0002127014
## n=20, result=0.7100538, error=0.0001734994
## n=21, result=0.7100844, error=0.0001429188
## n=22, result=0.7101085, error=0.0001187833
## n=23, result=0.7101278, error=9.953088e-05Jadi, dengan metode Trapezoida dibutuhkan n=23.
Metode Simpson
simpson_n <- function(ftn, a, b, n = 100) {
n <- max(c(2*(n %/% 2), 4))
h <- (b-a)/n
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))
return(S)
}f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
while(err>tol){
res_simp <- simpson_n(f,0,1,n = n)
err <- abs(res_simp-exact_value)
cat("n=",n,", result=",res_simp,", error=",err,"\n",sep = "")
n=n+1
if(n==1000){break}
}## n=4, result=0.7683974, error=0.05817006
## n=5, result=0.7683974, error=0.05817006
## n=6, result=0.7497082, error=0.03948094
## n=7, result=0.7497082, error=0.03948094
## n=8, result=0.728112, error=0.01788467
## n=9, result=0.728112, error=0.01788467
## n=10, result=0.7188088, error=0.008581532
## n=11, result=0.7188088, error=0.008581532
## n=12, result=0.7147289, error=0.004501593
## n=13, result=0.7147289, error=0.004501593
## n=14, result=0.712782, error=0.002554712
## n=15, result=0.712782, error=0.002554712
## n=16, result=0.711774, error=0.001546662
## n=17, result=0.711774, error=0.001546662
## n=18, result=0.7112144, error=0.0009870756
## n=19, result=0.7112144, error=0.0009870756
## n=20, result=0.7108852, error=0.0006578651
## n=21, result=0.7108852, error=0.0006578651
## n=22, result=0.7106819, error=0.0004545628
## n=23, result=0.7106819, error=0.0004545628
## n=24, result=0.7105511, error=0.0003237819
## n=25, result=0.7105511, error=0.0003237819
## n=26, result=0.710464, error=0.0002366803
## n=27, result=0.710464, error=0.0002366803
## n=28, result=0.7104042, error=0.0001769138
## n=29, result=0.7104042, error=0.0001769138
## n=30, result=0.7103621, error=0.0001348299
## n=31, result=0.7103621, error=0.0001348299
## n=32, result=0.7103318, error=0.0001045198
## n=33, result=0.7103318, error=0.0001045198
## n=34, result=0.7103096, error=8.225045e-05Jadi, dengan metode Simpson dibutuhkan n=34.
Metode Gauss Quadrature
f = function(x){15*x^5 - 30*x^10 + 15*x^15}
exact_value=0.7102273
tol <- 0.0001
err <- 1
n = 4
library(pracma)
while(err>tol){
gL <- gaussLegendre(n = n,a = 0,1)
Ci <- gL$w # koefisien
xi <- gL$x # gauss point
res_gl <- sum(Ci * f(xi))
err <- abs(res_gl-exact_value)
cat("n=",n,", result=",res_gl,", error=",err,"\n",sep = "")
n=n+1
if(n==1000){break}
}## n=4, result=0.6792999, error=0.03092738
## n=5, result=0.707375, error=0.002852279
## n=6, result=0.7101422, error=8.507726e-05Jadi, dengan Metode Gauss Quadrature dibutuhkan n=6.
Perbandingan
Berikut perbandingan banyaknya sub interval n untuk mendapatkan \(E(x)\) dengan toleransi 0.0001.
data.frame(Metode=c("Trapezoida","Simpson","Gauss Quadrature"),"Interval"=c(23,34,6))## Metode Interval
## 1 Trapezoida 23
## 2 Simpson 34
## 3 Gauss Quadrature 6Terlihat bahwa Gauss Quadrature adalah metode terbaik karena hanya membutuhkan 6 sub interval atau iterasinya paling sedikit.
Jawaban Soal 2
Pdf dari distirbusi Eksponensial adalah \[f(x|\lambda)=\lambda \ e^{-\lambda x}\] Karena \(\lambda=2\), maka fungsinya dapat ditulis:
\[f(x)=2\ e^{-2x}\] Untuk \(x=4\), CDF-nya dapat diperoleh sebagai berikut: \[F(4)=\int_{0}^{4} f(x)\ dx=\int_{0}^{4} 2\ e^{-2x}\ dx\] Dengan kata lain, akan dicari luas daerah di bawah kurva \(f(x)=2\ e^{-2x}\) dari x=0 sampai x=4.
library(DescTools)
f2 = function(x){ 2*exp(-2*x) }
curve(f2,xlim=c(0,5),ylab="y",main=expression(paste("f(x) = 2 ",e^{-2*x})))
abline(c(h=0,v=0))
Shade(f2,breaks=c(0,4),col="green")
grid()
Secara exact CDF dari \(f(x|\lambda)=\lambda \ e^{-\lambda x}\) adalah \(F(x)=1-e^{-\lambda x}\)
Jadi untuk \(\lambda=2\) dan \(x=4\), diperoleh \(F(x)=1-e^{-8}\)
cdf = 1-exp(-8)
cdf## [1] 0.9996645Dengan nilai exact \(cdf=0.9996645\), berikut dibandingkan hasil yang diperoleh melalui metode Trapezoida dan Simpson untuk n = 100, Gauss Quadrature dengan n = 7, dan Monte Carlo dengan m = 10000.
# Fungsi untuk Metode Monte Carlo
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)
}Menghitung CDF
cdf = 0.9996645
g = function(x){ 2*exp(-2*x) }
# Metode Trapezoida
h1 = trapezoid(g,a=0,b=4,n=100)
er1 = abs(cdf-h1)
# Metode Simpson
h2 = simpson_n(g,a=0,b=4,n=100)
er2 = abs(cdf-h2)
# Metode Gauss Quadrature
gL = gaussLegendre(n = 7,a = 0,4)
Ci = gL$w # koefisien
xi = gL$x # gauss point
h3 = sum(Ci * g(xi))
er3 = abs(h3-cdf)
# Metode Monte Carlo
set.seed(101)
m=10000
h4=mc_integral(g,a=0,b=4,m=m)
er4 = abs(h4-cdf)
Metode = c("Trapezoida","Simpson","Gauss Quadrature","Monte Carlo")
n = c(100,100,7,m)
Hasil = c(h1,h2,h3,h4)
Error = c(er1,er2,er3,er4)
data.frame(Metode,n,Hasil,Error)## Metode n Hasil Error
## 1 Trapezoida 100 1.0001976 5.331349e-04
## 2 Simpson 100 0.9996648 2.646781e-07
## 3 Gauss Quadrature 7 0.9996645 1.619311e-08
## 4 Monte Carlo 10000 0.9917179 7.946581e-03Terlihat bahwa Metode Gauss Quadrature memberikan hasil yang terbaik karena memeiliki error yang paling kecil.