Tugas 10 Praktikum STA561

Fungsi user-defined

Fungsi Trapezoidal

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)
}

Fungsi 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)
}

Fungsi Gaussian (Adaptive) Quadrature

gauss <- function (fx, n, a, b ){
  
#Legendre order = n
gL <- gaussLegendre(n = n, a = a,b)

#mendefinisikan koefisien dan gauss point

Ci <- gL$w # koefisien
xi <- gL$x # gauss point

#Menghitung integral 
I <- sum(Ci * fx(xi))

return (list(I=I,Ci=Ci,xi=xi))
}

Fungsi Monte Carlo Integration

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)
}

Soal 1

Diketahui pdf (probability density function) dari distribusi Kumaraswamy adalah sebagai berikut

\[f(x|a,b)=abx^{(a−1)}(1−x^{(a)})^{(b−1)};x∈(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.
Jika diketahui nilai Exact dari \(E(X)=\frac{bΓ(1+\frac1a)Γ(b)}{Γ(1+\frac1a+b)}\) metode mana yang paling mendekati untuk menghitung nilai \(E(X)\) dengan toleransi 0.0001? Gunakan R!

Jawab Soal 1 Manual (a)

Jika diketahui pdf dari distribusi Kusmarawamy dengan a = 5 , b = 3, maka didapatkan

\[f(x|5,3)=5⋅3x^{5−1}(1−x^5)^{3−1}=15x^4(1−x^5)^2\]

untuk \(E(X)\) yaitu

\(E(X)=\int_\infty^{-\infty}xf(x)dx\)

\(E(X)=\int_0^1x(15x^4(1−x^5)^2)dx\)

\(E(X)=\int_0^115x^5(1−x^5)^2dx\)

untuk selanjutnya dalam proses perhitungan, \(f(x)\) yang digunakan adalah fungsi \(E(x)\)

Metode Trapezoidal

\(\int_{a}^{b}f(x)dx \approx T_{n} = \frac{h}{2} [f(x_{0})+2f(x_{1})+2f(x_{2})+\ldots+2f(x_{n-1})+f(x_{n})]\)

Untuk n=4 maka

\(\int_{0}^{1} 15x^{5} (1-x^{5})^ {2}dx \approx T_{4} = \frac{h}{2} [f(x_{0})+2f(x_{1})+2f(x_{2})+2f(x_{3})+f(x_{4})]\)

Panjang antar sub interval \(h\) adalah

\(h=\frac{b-a}{n}=\frac{1-0}{4}=\frac{1}{4}= 0.25\)

\(f(x)\) pada setiap \(xi\) \((x1, x2, x3, x4)\)

\(f(x_{0})=f(0)=15 \sf x 0^{5} (1-0^{5})^ {2}=0\)

\(f(x_{1})=f(0.25)=15 \sf x 0.25^{5} (1-0.25^{5})^ {2}= 0.014619841\)

\(f(x_{3})=f(0.75)=15 {\sf x} 0.75^{5} (1-0.75^{5})^ {2}= 2.070616786\)

\(f(x_{4})=f(1)=15 {\sf x} 1^{5} (1-1^{5})^ {2}= 0\)

Maka \(T4\) dengan \(h=0.25\) adalah:

\(T_{4} = \frac{0.25}{2} [0 + (2 {\sf x} 0.014619841) + (2 {\sf x }0.439910889) + (2{\sf x} 2.070616786)+0]=0.631286879\)

Sehingga

\(\int_{0}^{1} 15x^{5} (1-x^{5})^ {2}dx \approx T_{4} = 0,63128\)

Metode Simpson

\(\int_{a}^{b}f(x)dx \approx S_{n} = \frac{h}{3} [f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\ldots+4f(x_{n-1})+f(x_{n})]\)

untuk n = 4

\(\int_{0}^{1}f(x)dx \approx S_{4} = \frac{h}{3} [f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+f(x_{4}))]\)

Panjang antar sub interval \(h\) adalah

\(h=\frac{b-a}{n}=\frac{1-0}{4}=\frac{1}{4}= 0.25\)

\(f(x)\) pada setiap \(xi\) \((x1, x2, x3, x4)\):

\(f(x_{0})=f(0)=15 \sf x 0^{5} (1-0^{5})^ {2}=0\)

\(f(x_{1})=f(0.25)=15 \sf x 0.25^{5} (1-0.25^{5})^ {2}= 0.014620\)

\(f(x_{2})=f(0.50)=15 {\sf x} 0.50^{5} (1-0.50^{5})^ {2}= 0.439911\)

\(f(x_{3})=f(0.75)=15 {\sf x} 0.75^{5} (1-0.75^{5})^ {2}= 2.070617\)

\(f(x_{4})=f(1)=15 {\sf x} 1^{5} (1-1^{5})^ {2}= 0\)

Maka \(S4\) dengan \(h=0.25\) adalah:

\(S_{4} = \frac{0.25}{2} [0 + (2 {\sf x} 0,014620) + (2 {\sf x } 0,439911) + (2{\sf x} 0,070617)+0]= 0.768397357\)

Metode Gauss Quadrature

Bentuk umum dari Gauss quadrature

\(\int_{a}^{b}f(x)dx \approx I= \sum_{i=1}^{n}C_{i}f(x_{i})\)

Untuk n=4 maka

\(\int_{a}^{b}f(x)dx \approx I= \sum_{i=1}^{4}C_{i}f(x_{i})\)

Mengubah domain fungsi integral ke [-1, 1]

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

dan

\(dx=\frac{1}{2}(b-a)dt = \frac{1}{2}dt\)

Sehingga

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

Menggunakan nilai koefisien dan titik gauss saat n=4

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

\(I= 0.3478548f(-0.8611363)+0.6521452f(-0.3399810)+0.6521452f(0.3399810)+0.3478548f(0.8611363)\)

Nilai \(f(x)\)

\(f(−0.8611363) = 0.0000042\)

\(f(−0,3399810) = 0.0189947\)

\(f(0,3399810) = 0.4940581\)

\(f(0.8611363) = 0.1662429\)

Maka nilai I adalah:

\(I = 0.3478548 \cdot 0.0000042+0.6521452 \cdot 0.0189947+0.6521452 \cdot 0.4940581 + 0.3478548 \cdot 0.1662429 = 0.679299903\)

Jawab Soal 1 dengan R (b)

Mendefinisikan fungsi E(x)

f <- function(x){
  15*x^(5)*((1-x^(5))^2)
}

Nilai Exact

Secara exact nilai E(X) adalah sebagai berikut:

a <- 5
b <- 3
(exact <- b*gamma(1+(1/a))*gamma(b)/(gamma(1+(1/a)+3)))

[1] 0.7102273

Atau dengan fungsi integrate ;

Fun <- function (fx, a, b) {
  integrate(fx, a,b)
}

Fun(f, 0,1)

0.7102273 with absolute error < 7.9e-15

\(E(X) = \frac{b\Gamma (1+ \frac{1}{a})\Gamma(b)}{\Gamma(1+ \frac{1}{a}+b)} = \frac{b\Gamma (1+ \frac{1}{5})\Gamma(3)}{\Gamma(1+ \frac{1}{5}+3)} = 0,7102273\)

Metode Trapezoidal

# Menggunakan user defined function
trapezoid(f,0,1,n = 4)

[1] 0.6312869

# Menggunakan fungsi trapzfun dari package pracma
library(pracma)
trapzfun(f,0,1, maxit = 4)

$value
[1] 0.7098069

$iter
[1] 4

$rel.err
[1] 0.005901199

Metode Simpson

simpson_n(f,0,1,n = 4)

[1] 0.7683974

Metode Gauss Quadrature

# Menggunakan R dengan fungsi gaussLegendre dengan domain [−1,1]
ft <- function(t){
  (15/2)*((1/2*(t+1))^5)*(1-(1/2*(t+1))^5)^2
}

domaingauss <- gauss(ft, 4, -1,1)
domaingauss

$I
[1] 0.6792999

$Ci
[1] 0.3478548 0.6521452 0.6521452 0.3478548

$xi
[1] -0.8611363 -0.3399810  0.3399810  0.8611363

gL <- gaussLegendre(n = 4,a = -1,1)

Ci <- gL$w # koefisien
Ci

[1] 0.3478548 0.6521452 0.6521452 0.3478548

xi <- gL$x # gauss point
xi

[1] -0.8611363 -0.3399810  0.3399810  0.8611363

I <- sum(Ci * ft(xi))
I

[1] 0.6792999

Soal 2

Diketahui pdf (probability density function) dari distirbusi Eksponensial adalah sebagai berikut

\(f(x|\lambda) = \lambda exp^{-\lambda x} ; x \in (0, \infty)\)

Jika diketahui \(λ = 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.

Formula untuk mencari cdf dari suatu distribusi adalah sebagai berikut :

\(F(x) = \int_{-\infty}^{x} f(t) dt\)

sehingga, cdf dari distribusi eksponensial:

\(F(x) = \int_{-\infty}^{x} \lambda exp^{-\lambda t} dt = F(x|\lambda)=1 - e^{-\lambda x}\)

Nilai Exact

Mendefinisikan fungsi f(x) untuk \(λ = 2\)

\(f(x|2) = 2 exp^{-2 x}\)

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

Melalui fungsi integrate yang telah dibuat

Fx_eksponensial <- Fun(fe,0,4)
Fx_eksponensial

0.9996645 with absolute error < 1.1e-14

Atau nilai tersebut juga dapat diperoleh dengan menggunakan CDF eksponensial

Fekspo <- 1 - exp(-2*4)
Fekspo

[1] 0.9996645

Selanjutnya akan dilakukan perhitungan dengan menggunakan metode integral numerik dengan menetapkan n = 100, m = 100 dan toleransi = 0,0001

Metode Trapezoidal

trapezoid(fe,0,4,n=100)

[1] 1.000198

#Menggunakan fungsi trapzfun dari package pracma
library(pracma)
trapzfun(fe,0,4, maxit = 100)

$value
[1] 0.9996646

$iter
[1] 14

$rel.err
[1] 5.958465e-08

Metode Trapezoidal

simpson_n(fe,0,4,n = 100)

[1] 0.9996648

#Menghitung hingga hasil integral mendekati hasil exact-nya dengan toleransi 0.0001

exact_value= 0.9996645
tol <-  0.0001
err <- 1
n = 100

while(err>tol){
  res_simp <- simpson_n(fe,0,4,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=100, result=0.9996648, error=2.646781e-07

Metode Gauss Quadrature

Mengubah domain fungsi integral ke [-1, 1]

\(\int_{0}^{4} 2exp ^{-2 x } dx = \int_{-1}^{1} 2exp ^{-2 (2[t+1])} 2dt = \int_{-1}^{1} 4exp ^{(-4t-4)}dt=\int_{-1}^{1} 4exp ^{(-4(t+1))}dt\)

ft <- function(t){
  4*exp(-4*(t+1))
}

gL <- gaussLegendre(n = 100, a = -1,1)

Ci <- gL$w # koefisien
Ci

 [1] 0.0007346345 0.0017093927 0.0026839254 0.0036559612 0.0046244501
 [6] 0.0055884280 0.0065469485 0.0074990733 0.0084438715 0.0093804197
[11] 0.0103078026 0.0112251140 0.0121314577 0.0130259479 0.0139077107
[16] 0.0147758845 0.0156296211 0.0164680862 0.0172904606 0.0180959407
[21] 0.0188837396 0.0196530875 0.0204032326 0.0211334421 0.0218430024
[26] 0.0225312203 0.0231974232 0.0238409603 0.0244612027 0.0250575445
[31] 0.0256294029 0.0261762192 0.0266974592 0.0271926134 0.0276611982
[36] 0.0281027557 0.0285168543 0.0289030896 0.0292610841 0.0295904881
[41] 0.0298909796 0.0301622651 0.0304040795 0.0306161866 0.0307983790
[46] 0.0309504789 0.0310723374 0.0311638357 0.0312248843 0.0312554235
[51] 0.0312554235 0.0312248843 0.0311638357 0.0310723374 0.0309504789
[56] 0.0307983790 0.0306161866 0.0304040795 0.0301622651 0.0298909796
[61] 0.0295904881 0.0292610841 0.0289030896 0.0285168543 0.0281027557
[66] 0.0276611982 0.0271926134 0.0266974592 0.0261762192 0.0256294029
[71] 0.0250575445 0.0244612027 0.0238409603 0.0231974232 0.0225312203
 [ reached getOption("max.print") -- omitted 25 entries ]

xi <- gL$x # gauss point
xi

 [1] -0.99971373 -0.99849195 -0.99629513 -0.99312494 -0.98898440 -0.98387754
 [7] -0.97780936 -0.97078578 -0.96281365 -0.95390078 -0.94405587 -0.93328854
[13] -0.92160930 -0.90902957 -0.89556164 -0.88121868 -0.86601469 -0.84996453
[19] -0.83308388 -0.81538924 -0.79689789 -0.77762791 -0.75759812 -0.73682809
[25] -0.71533812 -0.69314920 -0.67028302 -0.64676191 -0.62260886 -0.59784747
[31] -0.57250193 -0.54659701 -0.52015802 -0.49321079 -0.46578165 -0.43789740
[37] -0.40958529 -0.38087298 -0.35178853 -0.32236034 -0.29261719 -0.26258812
[43] -0.23230248 -0.20178986 -0.17108008 -0.14020314 -0.10918920 -0.07806858
[49] -0.04687168 -0.01562898  0.01562898  0.04687168  0.07806858  0.10918920
[55]  0.14020314  0.17108008  0.20178986  0.23230248  0.26258812  0.29261719
[61]  0.32236034  0.35178853  0.38087298  0.40958529  0.43789740  0.46578165
[67]  0.49321079  0.52015802  0.54659701  0.57250193  0.59784747  0.62260886
[73]  0.64676191  0.67028302  0.69314920
 [ reached getOption("max.print") -- omitted 25 entries ]

I <- sum(Ci * ft(xi))
I

[1] 0.9996645

Metode Monte Carlo

set.seed(123)
mc_integral(fe,a=0,b=4, m=100)

[1] 0.9330439

Berdasarkan keempat metode yang dilakukan, metode gauss quadrature memiliki nilai error terkecil jika dibandingkan metode lainnya. Sehingga dalam hal ini, metode gauss quadrature merupakan metode terbaik dalam melakukan perhitungan integral numerik. Sebelumnya juga dilakukan percobaan pada n yang berbeda namun tidak ditampilkan, dari beberapa percobaan yang dilakukan, semakin besar n yang digunakan, akan memberikan presisi nilai integral yang semakin baik.