Integral Numerik

Soal 1

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

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

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

Poin (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)

Jawab:

A. Metode Trapezoidal (n=4)

Berikut merupakan pdf (probability density function) dari distirbusi Kumaraswamy :

Untuk menghitung nilai harapan atau \(E(X)\) digunakan rumus:

Selanjutnya, dihitung panjang setiap sub interval sebagai berikut:

Karena n=4, maka bentuk dari Metode Trapezoidal adalah

Selanjutnya, dihitung nilai dari fungsi \(f(x)\) untuk setiap \(x_i\):

Berdasarkan hasil di atas diperoleh nilai \(h=0.25\), sehingga:

B. Metode Simpson (n=4)

Berikut merupakan pdf (probability density function) dari distirbusi Kumaraswamy :

Untuk menghitung nilai harapan atau \(E(X)\) digunakan rumus:

Selanjutnya, dihitung panjang setiap sub interval sebagai berikut:

Karena n=4, maka bentuk dari Metode Simpson adalah

Selanjutnya, dihitung nilai dari fungsi \(f(x)\) untuk setiap \(x_i\):

Berdasarkan hasil di atas diperoleh nilai \(h=0.25\), sehingga:

C. Metode Four-Point Gaussian Quadrature

Berikut merupakan pdf (probability density function) dari distirbusi Kumaraswamy :

Untuk menghitung nilai harapan atau \(E(X)\) digunakan rumus:

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:

Koefisien \(C_i\) merupakan bobot dan \(x_i\) merupakan titik Gauss yang berada dalam interval [1,1]. Oleh karena itu, dilakukan transformasi sehingga domain integral nya menjadi [1,1]:

dan

Sehingga, rumus untuk menghitung nilai harapan atau \(E(X)\) menjadi:

Selanjutnya, dihitung nilai dari fungsi \(f(t)\) untuk setiap \(t_i\), dengan \(t_i\) adalah titik gauss:

Berdasarkan hasil di atas, di peroleh nilai \(E(X)\) dengan menggunakan metode four-point gauss quadrature :

Hasil Integral 3 Metode

Metode <- c("Trapezoidal (n=4)",
            "Simpson (n=4)",
            "Four-Point Gaussian Quadrature")

Hasil_Integral<-c(0.631286879,
                  0.7683973572,
                  0.679299932)

Tabel_perbandingan <- data.frame(Metode,Hasil_Integral)
datatable(Tabel_perbandingan)

Poin (b)

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

Jawab:

Menghitung nilai eksak:

Metode Trapezoidal

library(pracma)

#user-defined function
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)*(1-(x^5))^2
}

exact_value=0.7102272727
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.07894039
## n=5, result=0.6738123, error=0.03641494
## n=6, result=0.6914928, error=0.01873451
## n=7, result=0.6997126, error=0.0105147
## n=8, result=0.7039057, error=0.006321576
## n=9, result=0.7062116, error=0.004015628
## n=10, result=0.7075597, error=0.002667566
## n=11, result=0.7083885, error=0.001838794
## n=12, result=0.7089199, error=0.001307412
## n=13, result=0.7092729, error=0.0009543875
## n=14, result=0.7095147, error=0.0007126217
## n=15, result=0.7096846, error=0.0005426466
## n=16, result=0.7098069, error=0.0004203769
## n=17, result=0.7098966, error=0.0003306343
## n=18, result=0.7099637, error=0.0002635799
## n=19, result=0.7100146, error=0.0002126741
## n=20, result=0.7100538, error=0.0001734721
## n=21, result=0.7100844, error=0.0001428915
## n=22, result=0.7101085, error=0.000118756
## n=23, result=0.7101278, error=9.950358e-05

Jadi, dengan menggunakan Metode Trapezoidal diperoleh nilai \(E(X)=0.7101278\) dengan nilai \(error=9.950358e-05\) dan \(n=23\).

Metode Simpson

#user-defined function
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)*(1-(x^5))^2
}

exact_value=0.7102272727
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==300){
    break
  }
  
}
## n=4, result=0.7683974, error=0.05817008
## n=5, result=0.7683974, error=0.05817008
## n=6, result=0.7497082, error=0.03948097
## n=7, result=0.7497082, error=0.03948097
## n=8, result=0.728112, error=0.0178847
## n=9, result=0.728112, error=0.0178847
## n=10, result=0.7188088, error=0.008581559
## n=11, result=0.7188088, error=0.008581559
## n=12, result=0.7147289, error=0.00450162
## n=13, result=0.7147289, error=0.00450162
## n=14, result=0.712782, error=0.002554739
## n=15, result=0.712782, error=0.002554739
## n=16, result=0.711774, error=0.001546689
## n=17, result=0.711774, error=0.001546689
## n=18, result=0.7112144, error=0.0009871029
## n=19, result=0.7112144, error=0.0009871029
## n=20, result=0.7108852, error=0.0006578924
## n=21, result=0.7108852, error=0.0006578924
## n=22, result=0.7106819, error=0.0004545901
## n=23, result=0.7106819, error=0.0004545901
## n=24, result=0.7105511, error=0.0003238092
## n=25, result=0.7105511, error=0.0003238092
## n=26, result=0.710464, error=0.0002367076
## n=27, result=0.710464, error=0.0002367076
## n=28, result=0.7104042, error=0.0001769411
## n=29, result=0.7104042, error=0.0001769411
## n=30, result=0.7103621, error=0.0001348572
## n=31, result=0.7103621, error=0.0001348572
## n=32, result=0.7103318, error=0.0001045471
## n=33, result=0.7103318, error=0.0001045471
## n=34, result=0.7103096, error=8.227775e-05

Jadi, dengan menggunakan Metode Simpson diperoleh nilai \(E(X)=0.7103096\) dengan nilai \(error=8.227775e-05\) dan \(n=34\).

Metode Gauss Quadrature

#integral 
ga<-function(t){
  (15/2)*((0.5*(t+1))^5)*(1-((0.5*(t+1))^5))^2
}

exact_value=0.7102272727
tol <-  0.0001
err <- 1
n = 4

while(err>tol){
  
  gL <- gaussLegendre(n = n, a = -1,1)

 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  
  res_gl <- sum(Ci * ga(xi))

  err <- abs(res_gl-exact_value)
  
  cat("n=",n,", result=",res_gl,", error=",err,"\n",sep = "")
  
  n=n+1
  if(n==300){
    break
  }
  
}
## n=4, result=0.6792999, error=0.03092735
## n=5, result=0.707375, error=0.002852251
## n=6, result=0.7101422, error=8.504996e-05

Jadi, dengan menggunakan Metode Gauss Quadrature diperoleh nilai \(E(X)=0.7101422\) dengan nilai \(error=8.504996e-05\) dan \(n=6\).

Perbandingan 3 Metode

Berikut merupakan perbandingan hasil perhitungan menggunakan 3 metode di atas dengan pendekatan nilai eksak \(E(X)=0.7102272727\) dan \(error=0.0001\), diperoleh nilai E(X) untuk masing-masing metode sebagai berikut:

Metode <- c("Trapezoidal",
            "Simpson",
            "Gauss Quadrature")

Hasil_Integral<-c(0.7101278,
                  0.7103096,
                  0.7101422)

Nilai_Error<-c(9.950358e-05,
               8.227775e-05,
               8.504996e-05)

n<-c(23,34,6)

Tabel_perbandingan1 <- data.frame(Metode,Hasil_Integral,Nilai_Error,n)
datatable(Tabel_perbandingan1)

Berdasarkan nilai error dan \(E(x)\) dari masing-masing metode, dapat dilihat bahwa metode yang menghasilkan nilai paling mendekati nilai eksak dari \(E(X)\) adalah Metode Simpson, karena metode ini menghasilkan nilai error paling kecil dibandingkan dengan metode lainnya.

Jadi, dapat disimpulkan bahwa metode integral numerik terbaik untuk menghitung nilai \(E(X)\) adalah Metode Simpson.

Soal 2

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

\(f(x|a,b)=\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.

Jawab:

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

CDF dari suatu distribusi didefinisikan sebagai berikut:

Karena nilai \(x\) dari distribusi eksponensial harus lebih besar dari nol, maka CDF dari distribusi exponential(2) untuk \(x=4\) dapat dituliskan sebagai

Membuat fungsi di R berdasarkan hasil di atas:

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

1. Metode Trapezoidal

Berikut merupakan user-defined function untuk Metode Trapezoidal:

#user-defined function
trapezoid1 <- function(ftn, a, b, n = 200) {
     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)
}

Dengan menggunakan user-defined function Metode Trapezoidal untuk \(n=200\), diperoleh CDF dari Eksponensial (2) dengan \(x=4\) :

hasil1<-trapezoid1(g,0,4,n = 200)
cat("nilai CDF=", hasil1)
## nilai CDF= 0.9997978

Selanjutnya, dilakukan perhitungan CDF sehingga hasilnya mendekati nilai eksak berikut:

nilai.eksak<-integrate(g, lower = 0, upper = 4)
nilai.eksak
## 0.9996645 with absolute error < 1.1e-14

dengan menggunakan batas error 0.00001, diperoleh:

exact_value=0.9996645
tol <-  0.00001
err <- 1
n = 10

while(err>tol){
  res_trap <- trapezoid(g,0,4,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=10, result=1.05242, error=0.05275531
## n=11, result=1.043343, error=0.04367879
## n=12, result=1.036418, error=0.03675326
## n=13, result=1.031015, error=0.0313503
## n=14, result=1.026719, error=0.0270549
## n=15, result=1.023249, error=0.02358421
## n=16, result=1.020405, error=0.02074012
## n=17, result=1.018045, error=0.01838055
## n=18, result=1.016066, error=0.0164015
## n=19, result=1.01439, error=0.0147254
## n=20, result=1.012958, error=0.01329349
## n=21, result=1.011725, error=0.01206057
## n=22, result=1.010656, error=0.01099142
## n=23, result=1.009723, error=0.01005831
## n=24, result=1.008904, error=0.009239095
## n=25, result=1.00818, error=0.008515985
## n=26, result=1.007539, error=0.00787452
## n=27, result=1.006967, error=0.007302862
## n=28, result=1.006456, error=0.006791242
## n=29, result=1.005996, error=0.006331541
## n=30, result=1.005581, error=0.005916966
## n=31, result=1.005206, error=0.005541802
## n=32, result=1.004866, error=0.005201208
## n=33, result=1.004556, error=0.004891065
## n=34, result=1.004272, error=0.00460785
## n=35, result=1.004013, error=0.004348533
## n=36, result=1.003775, error=0.004110501
## n=37, result=1.003556, error=0.003891487
## n=38, result=1.003354, error=0.003689518
## n=39, result=1.003167, error=0.00350287
## n=40, result=1.002995, error=0.003330033
## n=41, result=1.002834, error=0.003169677
## n=42, result=1.002685, error=0.003020629
## n=43, result=1.002546, error=0.00288185
## n=44, result=1.002417, error=0.002752418
## n=45, result=1.002296, error=0.002631513
## n=46, result=1.002183, error=0.002518402
## n=47, result=1.002077, error=0.002412428
## n=48, result=1.001978, error=0.002313005
## n=49, result=1.001884, error=0.002219603
## n=50, result=1.001796, error=0.002131746
## n=51, result=1.001714, error=0.002049003
## n=52, result=1.001635, error=0.001970985
## n=53, result=1.001562, error=0.001897339
## n=54, result=1.001492, error=0.001827745
## n=55, result=1.001426, error=0.00176191
## n=56, result=1.001364, error=0.001699569
## n=57, result=1.001305, error=0.001640479
## n=58, result=1.001249, error=0.001584417
## n=59, result=1.001196, error=0.001531181
## n=60, result=1.001145, error=0.001480583
## n=61, result=1.001097, error=0.001432452
## n=62, result=1.001051, error=0.001386631
## n=63, result=1.001007, error=0.001342973
## n=64, result=1.000966, error=0.001301345
## n=65, result=1.000926, error=0.001261623
## n=66, result=1.000888, error=0.001223692
## n=67, result=1.000852, error=0.001187446
## n=68, result=1.000817, error=0.001152787
## n=69, result=1.000784, error=0.001119624
## n=70, result=1.000752, error=0.001087871
## n=71, result=1.000722, error=0.00105745
## n=72, result=1.000693, error=0.001028287
## n=73, result=1.000665, error=0.001000315
## n=74, result=1.000638, error=0.000973468
## n=75, result=1.000612, error=0.0009476878
## n=76, result=1.000587, error=0.0009229182
## n=77, result=1.000564, error=0.0008991072
## n=78, result=1.000541, error=0.000876206
## n=79, result=1.000519, error=0.0008541686
## n=80, result=1.000497, error=0.0008329523
## n=81, result=1.000477, error=0.0008125168
## n=82, result=1.000457, error=0.0007928242
## n=83, result=1.000438, error=0.0007738389
## n=84, result=1.00042, error=0.0007555275
## n=85, result=1.000402, error=0.0007378584
## n=86, result=1.000385, error=0.000720802
## n=87, result=1.000369, error=0.0007043303
## n=88, result=1.000353, error=0.0006884168
## n=89, result=1.000338, error=0.0006730365
## n=90, result=1.000323, error=0.000658166
## n=91, result=1.000308, error=0.000643783
## n=92, result=1.000294, error=0.0006298663
## n=93, result=1.000281, error=0.0006163961
## n=94, result=1.000268, error=0.0006033534
## n=95, result=1.000255, error=0.0005907204
## n=96, result=1.000243, error=0.00057848
## n=97, result=1.000231, error=0.0005666162
## n=98, result=1.00022, error=0.0005551136
## n=99, result=1.000208, error=0.0005439578
## n=100, result=1.000198, error=0.0005331349
## n=101, result=1.000187, error=0.0005226319
## n=102, result=1.000177, error=0.0005124362
## n=103, result=1.000167, error=0.000502536
## n=104, result=1.000157, error=0.00049292
## n=105, result=1.000148, error=0.0004835774
## n=106, result=1.000139, error=0.0004744979
## n=107, result=1.00013, error=0.0004656717
## n=108, result=1.000122, error=0.0004570896
## n=109, result=1.000113, error=0.0004487425
## n=110, result=1.000105, error=0.000440622
## n=111, result=1.000097, error=0.00043272
## n=112, result=1.00009, error=0.0004250287
## n=113, result=1.000082, error=0.0004175406
## n=114, result=1.000075, error=0.0004102487
## n=115, result=1.000068, error=0.0004031462
## n=116, result=1.000061, error=0.0003962266
## n=117, result=1.000054, error=0.0003894836
## n=118, result=1.000047, error=0.0003829112
## n=119, result=1.000041, error=0.0003765039
## n=120, result=1.000035, error=0.0003702561
## n=121, result=1.000029, error=0.0003641625
## n=122, result=1.000023, error=0.0003582181
## n=123, result=1.000017, error=0.0003524181
## n=124, result=1.000011, error=0.0003467579
## n=125, result=1.000006, error=0.0003412329
## n=126, result=1, error=0.000335839
## n=127, result=0.9999951, error=0.0003305719
## n=128, result=0.9999899, error=0.0003254278
## n=129, result=0.9999849, error=0.0003204029
## n=130, result=0.99998, error=0.0003154935
## n=131, result=0.9999752, error=0.000310696
## n=132, result=0.9999705, error=0.0003060072
## n=133, result=0.9999659, error=0.0003014237
## n=134, result=0.9999614, error=0.0002969424
## n=135, result=0.9999571, error=0.0002925604
## n=136, result=0.9999528, error=0.0002882747
## n=137, result=0.9999486, error=0.0002840824
## n=138, result=0.9999445, error=0.000279981
## n=139, result=0.9999405, error=0.0002759677
## n=140, result=0.9999365, error=0.0002720401
## n=141, result=0.9999327, error=0.0002681958
## n=142, result=0.9999289, error=0.0002644324
## n=143, result=0.9999252, error=0.0002607477
## n=144, result=0.9999216, error=0.0002571395
## n=145, result=0.9999181, error=0.0002536057
## n=146, result=0.9999146, error=0.0002501442
## n=147, result=0.9999113, error=0.0002467531
## n=148, result=0.9999079, error=0.0002434306
## n=149, result=0.9999047, error=0.0002401747
## n=150, result=0.9999015, error=0.0002369837
## n=151, result=0.9998984, error=0.0002338558
## n=152, result=0.9998953, error=0.0002307895
## n=153, result=0.9998923, error=0.0002277832
## n=154, result=0.9998893, error=0.0002248352
## n=155, result=0.9998864, error=0.000221944
## n=156, result=0.9998836, error=0.0002191083
## n=157, result=0.9998808, error=0.0002163266
## n=158, result=0.9998781, error=0.0002135976
## n=159, result=0.9998754, error=0.0002109198
## n=160, result=0.9998728, error=0.0002082921
## n=161, result=0.9998702, error=0.0002057133
## n=162, result=0.9998677, error=0.000203182
## n=163, result=0.9998652, error=0.0002006972
## n=164, result=0.9998628, error=0.0001982577
## n=165, result=0.9998604, error=0.0001958624
## n=166, result=0.999858, error=0.0001935102
## n=167, result=0.9998557, error=0.0001912002
## n=168, result=0.9998534, error=0.0001889313
## n=169, result=0.9998512, error=0.0001867026
## n=170, result=0.999849, error=0.0001845131
## n=171, result=0.9998469, error=0.0001823618
## n=172, result=0.9998447, error=0.000180248
## n=173, result=0.9998427, error=0.0001781708
## n=174, result=0.9998406, error=0.0001761292
## n=175, result=0.9998386, error=0.0001741225
## n=176, result=0.9998366, error=0.00017215
## n=177, result=0.9998347, error=0.0001702108
## n=178, result=0.9998328, error=0.0001683042
## n=179, result=0.9998309, error=0.0001664294
## n=180, result=0.9998291, error=0.0001645858
## n=181, result=0.9998273, error=0.0001627727
## n=182, result=0.9998255, error=0.0001609893
## n=183, result=0.9998237, error=0.0001592352
## n=184, result=0.999822, error=0.0001575095
## n=185, result=0.9998203, error=0.0001558117
## n=186, result=0.9998186, error=0.0001541413
## n=187, result=0.999817, error=0.0001524976
## n=188, result=0.9998154, error=0.00015088
## n=189, result=0.9998138, error=0.0001492881
## n=190, result=0.9998122, error=0.0001477212
## n=191, result=0.9998107, error=0.0001461789
## n=192, result=0.9998092, error=0.0001446606
## n=193, result=0.9998077, error=0.0001431658
## n=194, result=0.9998062, error=0.0001416941
## n=195, result=0.9998047, error=0.000140245
## n=196, result=0.9998033, error=0.000138818
## n=197, result=0.9998019, error=0.0001374127
## n=198, result=0.9998005, error=0.0001360286
## n=199, result=0.9997992, error=0.0001346653
## n=200, result=0.9997978, error=0.0001333224
## n=201, result=0.9997965, error=0.0001319995
## n=202, result=0.9997952, error=0.0001306962
## n=203, result=0.9997939, error=0.0001294122
## n=204, result=0.9997926, error=0.0001281469
## n=205, result=0.9997914, error=0.0001269002
## n=206, result=0.9997902, error=0.0001256715
## n=207, result=0.999789, error=0.0001244606
## n=208, result=0.9997878, error=0.0001232671
## n=209, result=0.9997866, error=0.0001220908
## n=210, result=0.9997854, error=0.0001209311
## n=211, result=0.9997843, error=0.000119788
## n=212, result=0.9997832, error=0.0001186609
## n=213, result=0.999782, error=0.0001175497
## n=214, result=0.999781, error=0.0001164541
## n=215, result=0.9997799, error=0.0001153737
## n=216, result=0.9997788, error=0.0001143083
## n=217, result=0.9997778, error=0.0001132575
## n=218, result=0.9997767, error=0.0001122212
## n=219, result=0.9997757, error=0.0001111991
## n=220, result=0.9997747, error=0.0001101908
## n=221, result=0.9997737, error=0.0001091962
## n=222, result=0.9997727, error=0.0001082151
## n=223, result=0.9997717, error=0.000107247
## n=224, result=0.9997708, error=0.000106292
## n=225, result=0.9997698, error=0.0001053496
## n=226, result=0.9997689, error=0.0001044197
## n=227, result=0.999768, error=0.0001035021
## n=228, result=0.9997671, error=0.0001025965
## n=229, result=0.9997662, error=0.0001017028
## n=230, result=0.9997653, error=0.0001008207
## n=231, result=0.9997644, error=9.995e-05
## n=232, result=0.9997636, error=9.909056e-05
## n=233, result=0.9997627, error=9.824215e-05
## n=234, result=0.9997619, error=9.740461e-05
## n=235, result=0.9997611, error=9.657773e-05
## n=236, result=0.9997603, error=9.576134e-05
## n=237, result=0.9997595, error=9.495526e-05
## n=238, result=0.9997587, error=9.415932e-05
## n=239, result=0.9997579, error=9.337335e-05
## n=240, result=0.9997571, error=9.259719e-05
## n=241, result=0.9997563, error=9.183067e-05
## n=242, result=0.9997556, error=9.107362e-05
## n=243, result=0.9997548, error=9.032591e-05
## n=244, result=0.9997541, error=8.958737e-05
## n=245, result=0.9997534, error=8.885785e-05
## n=246, result=0.9997526, error=8.813722e-05
## n=247, result=0.9997519, error=8.742531e-05
## n=248, result=0.9997512, error=8.672201e-05
## n=249, result=0.9997505, error=8.602715e-05
## n=250, result=0.9997498, error=8.534062e-05
## n=251, result=0.9997492, error=8.466228e-05
## n=252, result=0.9997485, error=8.3992e-05
## n=253, result=0.9997478, error=8.332965e-05
## n=254, result=0.9997472, error=8.267511e-05
## n=255, result=0.9997465, error=8.202825e-05
## n=256, result=0.9997459, error=8.138896e-05
## n=257, result=0.9997453, error=8.075711e-05
## n=258, result=0.9997446, error=8.01326e-05
## n=259, result=0.999744, error=7.951531e-05
## n=260, result=0.9997434, error=7.890512e-05
## n=261, result=0.9997428, error=7.830194e-05
## n=262, result=0.9997422, error=7.770565e-05
## n=263, result=0.9997416, error=7.711615e-05
## n=264, result=0.999741, error=7.653333e-05
## n=265, result=0.9997405, error=7.59571e-05
## n=266, result=0.9997399, error=7.538736e-05
## n=267, result=0.9997393, error=7.482401e-05
## n=268, result=0.9997388, error=7.426695e-05
## n=269, result=0.9997382, error=7.371609e-05
## n=270, result=0.9997377, error=7.317134e-05
## n=271, result=0.9997371, error=7.263261e-05
## n=272, result=0.9997366, error=7.209981e-05
## n=273, result=0.9997361, error=7.157285e-05
## n=274, result=0.9997356, error=7.105166e-05
## n=275, result=0.999735, error=7.053614e-05
## n=276, result=0.9997345, error=7.002621e-05
## n=277, result=0.999734, error=6.952179e-05
## n=278, result=0.9997335, error=6.902281e-05
## n=279, result=0.999733, error=6.852918e-05
## n=280, result=0.9997325, error=6.804084e-05
## n=281, result=0.9997321, error=6.755769e-05
## n=282, result=0.9997316, error=6.707968e-05
## n=283, result=0.9997311, error=6.660673e-05
## n=284, result=0.9997306, error=6.613876e-05
## n=285, result=0.9997302, error=6.567571e-05
## n=286, result=0.9997297, error=6.521751e-05
## n=287, result=0.9997293, error=6.476409e-05
## n=288, result=0.9997288, error=6.431539e-05
## n=289, result=0.9997284, error=6.387133e-05
## n=290, result=0.9997279, error=6.343186e-05
## n=291, result=0.9997275, error=6.299691e-05
## n=292, result=0.9997271, error=6.256643e-05
## n=293, result=0.9997266, error=6.214034e-05
## n=294, result=0.9997262, error=6.17186e-05
## n=295, result=0.9997258, error=6.130113e-05
## n=296, result=0.9997254, error=6.088789e-05
## n=297, result=0.999725, error=6.047882e-05
## n=298, result=0.9997246, error=6.007386e-05
## n=299, result=0.9997242, error=5.967295e-05
## n=300, result=0.9997238, error=5.927605e-05
## n=301, result=0.9997234, error=5.88831e-05
## n=302, result=0.999723, error=5.849404e-05
## n=303, result=0.9997226, error=5.810883e-05
## n=304, result=0.9997222, error=5.772741e-05
## n=305, result=0.9997218, error=5.734974e-05
## n=306, result=0.9997215, error=5.697577e-05
## n=307, result=0.9997211, error=5.660544e-05
## n=308, result=0.9997207, error=5.623872e-05
## n=309, result=0.9997204, error=5.587555e-05
## n=310, result=0.99972, error=5.551588e-05
## n=311, result=0.9997197, error=5.515969e-05
## n=312, result=0.9997193, error=5.480691e-05
## n=313, result=0.999719, error=5.445751e-05
## n=314, result=0.9997186, error=5.411144e-05
## n=315, result=0.9997183, error=5.376866e-05
## n=316, result=0.9997179, error=5.342913e-05
## n=317, result=0.9997176, error=5.309281e-05
## n=318, result=0.9997173, error=5.275965e-05
## n=319, result=0.9997169, error=5.242963e-05
## n=320, result=0.9997166, error=5.210269e-05
## n=321, result=0.9997163, error=5.177881e-05
## n=322, result=0.999716, error=5.145793e-05
## n=323, result=0.9997156, error=5.114003e-05
## n=324, result=0.9997153, error=5.082508e-05
## n=325, result=0.999715, error=5.051302e-05
## n=326, result=0.9997147, error=5.020383e-05
## n=327, result=0.9997144, error=4.989748e-05
## n=328, result=0.9997141, error=4.959392e-05
## n=329, result=0.9997138, error=4.929312e-05
## n=330, result=0.9997135, error=4.899506e-05
## n=331, result=0.9997132, error=4.869969e-05
## n=332, result=0.9997129, error=4.840699e-05
## n=333, result=0.9997126, error=4.811692e-05
## n=334, result=0.9997123, error=4.782945e-05
## n=335, result=0.999712, error=4.754455e-05
## n=336, result=0.9997118, error=4.72622e-05
## n=337, result=0.9997115, error=4.698235e-05
## n=338, result=0.9997112, error=4.670498e-05
## n=339, result=0.9997109, error=4.643007e-05
## n=340, result=0.9997107, error=4.615757e-05
## n=341, result=0.9997104, error=4.588747e-05
## n=342, result=0.9997101, error=4.561974e-05
## n=343, result=0.9997099, error=4.535434e-05
## n=344, result=0.9997096, error=4.509125e-05
## n=345, result=0.9997093, error=4.483045e-05
## n=346, result=0.9997091, error=4.457191e-05
## n=347, result=0.9997088, error=4.43156e-05
## n=348, result=0.9997086, error=4.406149e-05
## n=349, result=0.9997083, error=4.380957e-05
## n=350, result=0.9997081, error=4.35598e-05
## n=351, result=0.9997078, error=4.331217e-05
## n=352, result=0.9997076, error=4.306664e-05
## n=353, result=0.9997073, error=4.282319e-05
## n=354, result=0.9997071, error=4.258181e-05
## n=355, result=0.9997068, error=4.234246e-05
## n=356, result=0.9997066, error=4.210513e-05
## n=357, result=0.9997064, error=4.186979e-05
## n=358, result=0.9997061, error=4.163642e-05
## n=359, result=0.9997059, error=4.140499e-05
## n=360, result=0.9997057, error=4.117549e-05
## n=361, result=0.9997054, error=4.09479e-05
## n=362, result=0.9997052, error=4.072219e-05
## n=363, result=0.999705, error=4.049834e-05
## n=364, result=0.9997048, error=4.027633e-05
## n=365, result=0.9997046, error=4.005615e-05
## n=366, result=0.9997043, error=3.983777e-05
## n=367, result=0.9997041, error=3.962117e-05
## n=368, result=0.9997039, error=3.940633e-05
## n=369, result=0.9997037, error=3.919324e-05
## n=370, result=0.9997035, error=3.898188e-05
## n=371, result=0.9997033, error=3.877222e-05
## n=372, result=0.9997031, error=3.856425e-05
## n=373, result=0.9997029, error=3.835795e-05
## n=374, result=0.9997027, error=3.81533e-05
## n=375, result=0.9997025, error=3.795029e-05
## n=376, result=0.9997022, error=3.774889e-05
## n=377, result=0.999702, error=3.75491e-05
## n=378, result=0.9997019, error=3.735089e-05
## n=379, result=0.9997017, error=3.715424e-05
## n=380, result=0.9997015, error=3.695915e-05
## n=381, result=0.9997013, error=3.676559e-05
## n=382, result=0.9997011, error=3.657355e-05
## n=383, result=0.9997009, error=3.638301e-05
## n=384, result=0.9997007, error=3.619396e-05
## n=385, result=0.9997005, error=3.600638e-05
## n=386, result=0.9997003, error=3.582025e-05
## n=387, result=0.9997001, error=3.563557e-05
## n=388, result=0.9997, error=3.545231e-05
## n=389, result=0.9996998, error=3.527046e-05
## n=390, result=0.9996996, error=3.509001e-05
## n=391, result=0.9996994, error=3.491095e-05
## n=392, result=0.9996992, error=3.473325e-05
## n=393, result=0.9996991, error=3.455691e-05
## n=394, result=0.9996989, error=3.43819e-05
## n=395, result=0.9996987, error=3.420823e-05
## n=396, result=0.9996985, error=3.403587e-05
## n=397, result=0.9996984, error=3.386481e-05
## n=398, result=0.9996982, error=3.369503e-05
## n=399, result=0.999698, error=3.352654e-05
## n=400, result=0.9996979, error=3.33593e-05
## n=401, result=0.9996977, error=3.319332e-05
## n=402, result=0.9996975, error=3.302857e-05
## n=403, result=0.9996974, error=3.286504e-05
## n=404, result=0.9996972, error=3.270273e-05
## n=405, result=0.999697, error=3.254162e-05
## n=406, result=0.9996969, error=3.23817e-05
## n=407, result=0.9996967, error=3.222296e-05
## n=408, result=0.9996966, error=3.206538e-05
## n=409, result=0.9996964, error=3.190895e-05
## n=410, result=0.9996963, error=3.175367e-05
## n=411, result=0.9996961, error=3.159953e-05
## n=412, result=0.9996959, error=3.14465e-05
## n=413, result=0.9996958, error=3.129458e-05
## n=414, result=0.9996956, error=3.114376e-05
## n=415, result=0.9996955, error=3.099403e-05
## n=416, result=0.9996953, error=3.084538e-05
## n=417, result=0.9996952, error=3.06978e-05
## n=418, result=0.9996951, error=3.055128e-05
## n=419, result=0.9996949, error=3.04058e-05
## n=420, result=0.9996948, error=3.026136e-05
## n=421, result=0.9996946, error=3.011795e-05
## n=422, result=0.9996945, error=2.997556e-05
## n=423, result=0.9996943, error=2.983418e-05
## n=424, result=0.9996942, error=2.969379e-05
## n=425, result=0.9996941, error=2.95544e-05
## n=426, result=0.9996939, error=2.941598e-05
## n=427, result=0.9996938, error=2.927854e-05
## n=428, result=0.9996936, error=2.914206e-05
## n=429, result=0.9996935, error=2.900653e-05
## n=430, result=0.9996934, error=2.887195e-05
## n=431, result=0.9996932, error=2.87383e-05
## n=432, result=0.9996931, error=2.860558e-05
## n=433, result=0.999693, error=2.847378e-05
## n=434, result=0.9996928, error=2.834289e-05
## n=435, result=0.9996927, error=2.82129e-05
## n=436, result=0.9996926, error=2.80838e-05
## n=437, result=0.9996925, error=2.795559e-05
## n=438, result=0.9996923, error=2.782826e-05
## n=439, result=0.9996922, error=2.770179e-05
## n=440, result=0.9996921, error=2.757619e-05
## n=441, result=0.999692, error=2.745144e-05
## n=442, result=0.9996918, error=2.732753e-05
## n=443, result=0.9996917, error=2.720447e-05
## n=444, result=0.9996916, error=2.708223e-05
## n=445, result=0.9996915, error=2.696082e-05
## n=446, result=0.9996913, error=2.684022e-05
## n=447, result=0.9996912, error=2.672043e-05
## n=448, result=0.9996911, error=2.660145e-05
## n=449, result=0.999691, error=2.648325e-05
## n=450, result=0.9996909, error=2.636585e-05
## n=451, result=0.9996907, error=2.624922e-05
## n=452, result=0.9996906, error=2.613337e-05
## n=453, result=0.9996905, error=2.601828e-05
## n=454, result=0.9996904, error=2.590396e-05
## n=455, result=0.9996903, error=2.579038e-05
## n=456, result=0.9996902, error=2.567755e-05
## n=457, result=0.9996901, error=2.556547e-05
## n=458, result=0.99969, error=2.545411e-05
## n=459, result=0.9996898, error=2.534349e-05
## n=460, result=0.9996897, error=2.523358e-05
## n=461, result=0.9996896, error=2.512439e-05
## n=462, result=0.9996895, error=2.50159e-05
## n=463, result=0.9996894, error=2.490812e-05
## n=464, result=0.9996893, error=2.480104e-05
## n=465, result=0.9996892, error=2.469464e-05
## n=466, result=0.9996891, error=2.458893e-05
## n=467, result=0.999689, error=2.44839e-05
## n=468, result=0.9996889, error=2.437954e-05
## n=469, result=0.9996888, error=2.427584e-05
## n=470, result=0.9996887, error=2.417281e-05
## n=471, result=0.9996886, error=2.407043e-05
## n=472, result=0.9996885, error=2.396871e-05
## n=473, result=0.9996884, error=2.386763e-05
## n=474, result=0.9996883, error=2.376718e-05
## n=475, result=0.9996882, error=2.366737e-05
## n=476, result=0.9996881, error=2.356819e-05
## n=477, result=0.999688, error=2.346963e-05
## n=478, result=0.9996879, error=2.337169e-05
## n=479, result=0.9996878, error=2.327437e-05
## n=480, result=0.9996877, error=2.317765e-05
## n=481, result=0.9996876, error=2.308153e-05
## n=482, result=0.9996875, error=2.298601e-05
## n=483, result=0.9996874, error=2.289108e-05
## n=484, result=0.9996873, error=2.279675e-05
## n=485, result=0.9996872, error=2.270299e-05
## n=486, result=0.9996871, error=2.260981e-05
## n=487, result=0.999687, error=2.251721e-05
## n=488, result=0.9996869, error=2.242517e-05
## n=489, result=0.9996868, error=2.23337e-05
## n=490, result=0.9996867, error=2.224279e-05
## n=491, result=0.9996867, error=2.215243e-05
## n=492, result=0.9996866, error=2.206262e-05
## n=493, result=0.9996865, error=2.197336e-05
## n=494, result=0.9996864, error=2.188464e-05
## n=495, result=0.9996863, error=2.179646e-05
## n=496, result=0.9996862, error=2.170881e-05
## n=497, result=0.9996861, error=2.162169e-05
## n=498, result=0.999686, error=2.153509e-05
## n=499, result=0.9996859, error=2.144902e-05
## n=500, result=0.9996859, error=2.136346e-05
## n=501, result=0.9996858, error=2.127841e-05
## n=502, result=0.9996857, error=2.119387e-05
## n=503, result=0.9996856, error=2.110983e-05
## n=504, result=0.9996855, error=2.102629e-05
## n=505, result=0.9996854, error=2.094325e-05
## n=506, result=0.9996854, error=2.08607e-05
## n=507, result=0.9996853, error=2.077864e-05
## n=508, result=0.9996852, error=2.069706e-05
## n=509, result=0.9996851, error=2.061596e-05
## n=510, result=0.999685, error=2.053534e-05
## n=511, result=0.999685, error=2.04552e-05
## n=512, result=0.9996849, error=2.037552e-05
## n=513, result=0.9996848, error=2.02963e-05
## n=514, result=0.9996847, error=2.021755e-05
## n=515, result=0.9996846, error=2.013926e-05
## n=516, result=0.9996846, error=2.006142e-05
## n=517, result=0.9996845, error=1.998403e-05
## n=518, result=0.9996844, error=1.990709e-05
## n=519, result=0.9996843, error=1.98306e-05
## n=520, result=0.9996843, error=1.975454e-05
## n=521, result=0.9996842, error=1.967893e-05
## n=522, result=0.9996841, error=1.960374e-05
## n=523, result=0.999684, error=1.952899e-05
## n=524, result=0.999684, error=1.945467e-05
## n=525, result=0.9996839, error=1.938077e-05
## n=526, result=0.9996838, error=1.930729e-05
## n=527, result=0.9996837, error=1.923423e-05
## n=528, result=0.9996837, error=1.916158e-05
## n=529, result=0.9996836, error=1.908935e-05
## n=530, result=0.9996835, error=1.901752e-05
## n=531, result=0.9996834, error=1.89461e-05
## n=532, result=0.9996834, error=1.887508e-05
## n=533, result=0.9996833, error=1.880446e-05
## n=534, result=0.9996832, error=1.873424e-05
## n=535, result=0.9996832, error=1.866441e-05
## n=536, result=0.9996831, error=1.859497e-05
## n=537, result=0.999683, error=1.852592e-05
## n=538, result=0.999683, error=1.845725e-05
## n=539, result=0.9996829, error=1.838897e-05
## n=540, result=0.9996828, error=1.832106e-05
## n=541, result=0.9996828, error=1.825353e-05
## n=542, result=0.9996827, error=1.818638e-05
## n=543, result=0.9996826, error=1.811959e-05
## n=544, result=0.9996826, error=1.805318e-05
## n=545, result=0.9996825, error=1.798712e-05
## n=546, result=0.9996824, error=1.792143e-05
## n=547, result=0.9996824, error=1.78561e-05
## n=548, result=0.9996823, error=1.779113e-05
## n=549, result=0.9996822, error=1.772651e-05
## n=550, result=0.9996822, error=1.766225e-05
## n=551, result=0.9996821, error=1.759833e-05
## n=552, result=0.999682, error=1.753476e-05
## n=553, result=0.999682, error=1.747154e-05
## n=554, result=0.9996819, error=1.740866e-05
## n=555, result=0.9996818, error=1.734612e-05
## n=556, result=0.9996818, error=1.728391e-05
## n=557, result=0.9996817, error=1.722204e-05
## n=558, result=0.9996817, error=1.71605e-05
## n=559, result=0.9996816, error=1.709929e-05
## n=560, result=0.9996815, error=1.703841e-05
## n=561, result=0.9996815, error=1.697786e-05
## n=562, result=0.9996814, error=1.691762e-05
## n=563, result=0.9996814, error=1.685771e-05
## n=564, result=0.9996813, error=1.679812e-05
## n=565, result=0.9996812, error=1.673884e-05
## n=566, result=0.9996812, error=1.667988e-05
## n=567, result=0.9996811, error=1.662123e-05
## n=568, result=0.9996811, error=1.656288e-05
## n=569, result=0.999681, error=1.650485e-05
## n=570, result=0.9996809, error=1.644712e-05
## n=571, result=0.9996809, error=1.638969e-05
## n=572, result=0.9996808, error=1.633257e-05
## n=573, result=0.9996808, error=1.627574e-05
## n=574, result=0.9996807, error=1.621921e-05
## n=575, result=0.9996807, error=1.616297e-05
## n=576, result=0.9996806, error=1.610703e-05
## n=577, result=0.9996806, error=1.605138e-05
## n=578, result=0.9996805, error=1.599601e-05
## n=579, result=0.9996804, error=1.594094e-05
## n=580, result=0.9996804, error=1.588615e-05
## n=581, result=0.9996803, error=1.583164e-05
## n=582, result=0.9996803, error=1.577741e-05
## n=583, result=0.9996802, error=1.572346e-05
## n=584, result=0.9996802, error=1.566978e-05
## n=585, result=0.9996801, error=1.561638e-05
## n=586, result=0.9996801, error=1.556326e-05
## n=587, result=0.99968, error=1.551041e-05
## n=588, result=0.99968, error=1.545782e-05
## n=589, result=0.9996799, error=1.54055e-05
## n=590, result=0.9996799, error=1.535345e-05
## n=591, result=0.9996798, error=1.530167e-05
## n=592, result=0.9996798, error=1.525014e-05
## n=593, result=0.9996797, error=1.519888e-05
## n=594, result=0.9996796, error=1.514787e-05
## n=595, result=0.9996796, error=1.509712e-05
## n=596, result=0.9996795, error=1.504663e-05
## n=597, result=0.9996795, error=1.499639e-05
## n=598, result=0.9996794, error=1.49464e-05
## n=599, result=0.9996794, error=1.489666e-05
## n=600, result=0.9996793, error=1.484717e-05
## n=601, result=0.9996793, error=1.479793e-05
## n=602, result=0.9996792, error=1.474893e-05
## n=603, result=0.9996792, error=1.470018e-05
## n=604, result=0.9996792, error=1.465167e-05
## n=605, result=0.9996791, error=1.46034e-05
## n=606, result=0.9996791, error=1.455536e-05
## n=607, result=0.999679, error=1.450757e-05
## n=608, result=0.999679, error=1.446001e-05
## n=609, result=0.9996789, error=1.441268e-05
## n=610, result=0.9996789, error=1.436559e-05
## n=611, result=0.9996788, error=1.431873e-05
## n=612, result=0.9996788, error=1.427209e-05
## n=613, result=0.9996787, error=1.422569e-05
## n=614, result=0.9996787, error=1.417951e-05
## n=615, result=0.9996786, error=1.413356e-05
## n=616, result=0.9996786, error=1.408783e-05
## n=617, result=0.9996785, error=1.404232e-05
## n=618, result=0.9996785, error=1.399703e-05
## n=619, result=0.9996785, error=1.395197e-05
## n=620, result=0.9996784, error=1.390712e-05
## n=621, result=0.9996784, error=1.386248e-05
## n=622, result=0.9996783, error=1.381806e-05
## n=623, result=0.9996783, error=1.377386e-05
## n=624, result=0.9996782, error=1.372987e-05
## n=625, result=0.9996782, error=1.368609e-05
## n=626, result=0.9996781, error=1.364252e-05
## n=627, result=0.9996781, error=1.359915e-05
## n=628, result=0.9996781, error=1.3556e-05
## n=629, result=0.999678, error=1.351305e-05
## n=630, result=0.999678, error=1.34703e-05
## n=631, result=0.9996779, error=1.342776e-05
## n=632, result=0.9996779, error=1.338542e-05
## n=633, result=0.9996778, error=1.334328e-05
## n=634, result=0.9996778, error=1.330134e-05
## n=635, result=0.9996778, error=1.325959e-05
## n=636, result=0.9996777, error=1.321805e-05
## n=637, result=0.9996777, error=1.31767e-05
## n=638, result=0.9996776, error=1.313554e-05
## n=639, result=0.9996776, error=1.309458e-05
## n=640, result=0.9996776, error=1.30538e-05
## n=641, result=0.9996775, error=1.301322e-05
## n=642, result=0.9996775, error=1.297283e-05
## n=643, result=0.9996774, error=1.293263e-05
## n=644, result=0.9996774, error=1.289261e-05
## n=645, result=0.9996774, error=1.285278e-05
## n=646, result=0.9996773, error=1.281314e-05
## n=647, result=0.9996773, error=1.277367e-05
## n=648, result=0.9996772, error=1.273439e-05
## n=649, result=0.9996772, error=1.26953e-05
## n=650, result=0.9996772, error=1.265638e-05
## n=651, result=0.9996771, error=1.261764e-05
## n=652, result=0.9996771, error=1.257908e-05
## n=653, result=0.999677, error=1.25407e-05
## n=654, result=0.999677, error=1.250249e-05
## n=655, result=0.999677, error=1.246446e-05
## n=656, result=0.9996769, error=1.24266e-05
## n=657, result=0.9996769, error=1.238891e-05
## n=658, result=0.9996769, error=1.23514e-05
## n=659, result=0.9996768, error=1.231406e-05
## n=660, result=0.9996768, error=1.227688e-05
## n=661, result=0.9996767, error=1.223988e-05
## n=662, result=0.9996767, error=1.220304e-05
## n=663, result=0.9996767, error=1.216637e-05
## n=664, result=0.9996766, error=1.212986e-05
## n=665, result=0.9996766, error=1.209352e-05
## n=666, result=0.9996766, error=1.205735e-05
## n=667, result=0.9996765, error=1.202133e-05
## n=668, result=0.9996765, error=1.198548e-05
## n=669, result=0.9996764, error=1.194979e-05
## n=670, result=0.9996764, error=1.191425e-05
## n=671, result=0.9996764, error=1.187888e-05
## n=672, result=0.9996763, error=1.184366e-05
## n=673, result=0.9996763, error=1.18086e-05
## n=674, result=0.9996763, error=1.17737e-05
## n=675, result=0.9996762, error=1.173895e-05
## n=676, result=0.9996762, error=1.170436e-05
## n=677, result=0.9996762, error=1.166992e-05
## n=678, result=0.9996761, error=1.163563e-05
## n=679, result=0.9996761, error=1.160149e-05
## n=680, result=0.9996761, error=1.15675e-05
## n=681, result=0.999676, error=1.153366e-05
## n=682, result=0.999676, error=1.149998e-05
## n=683, result=0.999676, error=1.146643e-05
## n=684, result=0.9996759, error=1.143304e-05
## n=685, result=0.9996759, error=1.139979e-05
## n=686, result=0.9996759, error=1.136669e-05
## n=687, result=0.9996758, error=1.133373e-05
## n=688, result=0.9996758, error=1.130092e-05
## n=689, result=0.9996758, error=1.126825e-05
## n=690, result=0.9996757, error=1.123572e-05
## n=691, result=0.9996757, error=1.120333e-05
## n=692, result=0.9996757, error=1.117108e-05
## n=693, result=0.9996756, error=1.113897e-05
## n=694, result=0.9996756, error=1.1107e-05
## n=695, result=0.9996756, error=1.107517e-05
## n=696, result=0.9996755, error=1.104348e-05
## n=697, result=0.9996755, error=1.101192e-05
## n=698, result=0.9996755, error=1.098049e-05
## n=699, result=0.9996754, error=1.094921e-05
## n=700, result=0.9996754, error=1.091805e-05
## n=701, result=0.9996754, error=1.088703e-05
## n=702, result=0.9996754, error=1.085614e-05
## n=703, result=0.9996753, error=1.082538e-05
## n=704, result=0.9996753, error=1.079476e-05
## n=705, result=0.9996753, error=1.076426e-05
## n=706, result=0.9996752, error=1.07339e-05
## n=707, result=0.9996752, error=1.070366e-05
## n=708, result=0.9996752, error=1.067355e-05
## n=709, result=0.9996751, error=1.064357e-05
## n=710, result=0.9996751, error=1.061371e-05
## n=711, result=0.9996751, error=1.058398e-05
## n=712, result=0.9996751, error=1.055438e-05
## n=713, result=0.999675, error=1.05249e-05
## n=714, result=0.999675, error=1.049554e-05
## n=715, result=0.999675, error=1.046631e-05
## n=716, result=0.9996749, error=1.04372e-05
## n=717, result=0.9996749, error=1.040821e-05
## n=718, result=0.9996749, error=1.037934e-05
## n=719, result=0.9996749, error=1.035059e-05
## n=720, result=0.9996748, error=1.032197e-05
## n=721, result=0.9996748, error=1.029346e-05
## n=722, result=0.9996748, error=1.026507e-05
## n=723, result=0.9996747, error=1.023679e-05
## n=724, result=0.9996747, error=1.020864e-05
## n=725, result=0.9996747, error=1.01806e-05
## n=726, result=0.9996747, error=1.015268e-05
## n=727, result=0.9996746, error=1.012487e-05
## n=728, result=0.9996746, error=1.009717e-05
## n=729, result=0.9996746, error=1.006959e-05
## n=730, result=0.9996745, error=1.004213e-05
## n=731, result=0.9996745, error=1.001477e-05
## n=732, result=0.9996745, error=9.98753e-06

Jadi, dengan menggunakan Metode Trapezoidal diperoleh nilai \(CDF=0.9996745\) dengan nilai \(error=9.98753e-06\) dan \(n=732\).

2. Metode Simpson

Berikut merupakan user-defined function untuk Metode Simpson:

#user-defined function
simpson_n1<- 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)
}

Dengan menggunakan user-defined function Metode Simpson untuk \(n=20\), diperoleh CDF dari Eksponensial (2) dengan \(x=4\) :

hasils1<-simpson_n1(g,0,4,n = 20)
cat("nilai CDF=", hasils1)
## nilai CDF= 0.999804

Selanjutnya, dilakukan perhitungan untuk nilai CDF sehingga hasilnya mendekati nilai eksak berikut:

nilai.eksak<-integrate(g, lower = 0, upper = 4)
nilai.eksak
## 0.9996645 with absolute error < 1.1e-14

dengan menggunakan batas error 0.00001, diperoleh nilai CDF:

exact_value=0.9996645
tol <- 0.00001
err <- 1
n = 5

while(err>tol){
  res_simp <- simpson_n(g,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=5, result=1.058815, error=0.05915075
## n=6, result=1.014089, error=0.01442435
## n=7, result=1.014089, error=0.01442435
## n=8, result=1.00462, error=0.004955555
## n=9, result=1.00462, error=0.004955555
## n=10, result=1.001777, error=0.002112444
## n=11, result=1.001777, error=0.002112444
## n=12, result=1.000706, error=0.00104161
## n=13, result=1.000706, error=0.00104161
## n=14, result=1.000234, error=0.0005699298
## n=15, result=1.000234, error=0.0005699298
## n=16, result=1.000002, error=0.000337077
## n=17, result=1.000002, error=0.000337077
## n=18, result=0.9998762, error=0.0002117416
## n=19, result=0.9998762, error=0.0002117416
## n=20, result=0.999804, error=0.0001395486
## n=21, result=0.999804, error=0.0001395486
## n=22, result=0.9997601, error=9.563677e-05
## n=23, result=0.9997601, error=9.563677e-05
## n=24, result=0.9997322, error=6.7705e-05
## n=25, result=0.9997322, error=6.7705e-05
## n=26, result=0.9997138, error=4.926109e-05
## n=27, result=0.9997138, error=4.926109e-05
## n=28, result=0.9997012, error=3.668996e-05
## n=29, result=0.9997012, error=3.668996e-05
## n=30, result=0.9996924, error=2.788524e-05
## n=31, result=0.9996924, error=2.788524e-05
## n=32, result=0.9996861, error=2.157112e-05
## n=33, result=0.9996861, error=2.157112e-05
## n=34, result=0.9996814, error=1.694843e-05
## n=35, result=0.9996814, error=1.694843e-05
## n=36, result=0.999678, error=1.350169e-05
## n=37, result=0.999678, error=1.350169e-05
## n=38, result=0.9996754, error=1.088964e-05
## n=39, result=0.9996754, error=1.088964e-05
## n=40, result=0.9996734, error=8.881142e-06

Jadi, dengan menggunakan Metode Simpson diperoleh nilai \(CDF=0.9996734\) dengan nilai \(error=8.881142e-06\) dan \(n=40\).

3. Metode Gauss Quadrature

Berdasarkan hasil di atas CDF untuk eksponensial (2) dengan x=4 dapat dihitung dengan menggunakan:

Karena koefisien Ci yang merupakan bobot dan xi merupakan titik Gauss berada dalam interval [???1,1] maka perlu dilakukan transformasi domain menjadi [-1,1]:

dan

Sehingga, integral untuk menghitung CDF menjadi:

Membuat fungsi di R berdasarkan hasil transformasi di atas:

g1<-function(x){
  4*exp(-4*(x+1))
}

Dengan menggunakan \(n=5\), diperoleh CDF dari Eksponensial (2) dengan \(x=4\) :

gL <- gaussLegendre(n = 5, a = -1,1)
Ci <- gL$w
xi <- gL$x 
I <- sum(Ci * g1(xi))
I
## [1] 0.9995785

Selanjutnya, dilakukan perhitungan untuk nilai CDF sehingga hasilnya mendekati nilai eksak berikut:

nilai.eksak<-integrate(g, lower = 0, upper = 4)
nilai.eksak
## 0.9996645 with absolute error < 1.1e-14

dengan menggunakan batas error 0.00001, diperoleh:

exact_value=0.9996645
tol <-  0.00001
err <- 1
n = 5

while(err>tol){
  
  gL <- gaussLegendre(n = n, a = -1,1)

 Ci <- gL$w # koefisien
 xi <- gL$x # gauss point
  
  res_gl <- sum(Ci * g1(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=5, result=0.9995785, error=8.59847e-05
## n=6, result=0.999662, error=2.471959e-06

Jadi, dengan menggunakan Metode Gauss Quadrature diperoleh nilai \(CDF=0.999662\) dengan nilai \(error=2.471959e-06\) dan \(n=6\).

4. Metode Monte Carlo

Berikut merupakan user-defined function untuk Metode Monte Carlo:

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

Dengan menggunakan \(m=11000\), diperoleh CDF dari Eksponensial (2) dengan \(x=4\) :

set.seed(123)
int_g1 <- mc_integral(g,a=0,b=4,m = 11000)
int_g1
## [1] 0.9952722

\[F(x=4|2)\approx0.9952722\]

Selanjutnya, dilakukan perhitungan untuk nilai CDF sehingga hasilnya mendekati nilai eksak berikut:

nilai.eksak<-integrate(g, lower = 0, upper = 4)
nilai.eksak
## 0.9996645 with absolute error < 1.1e-14

dengan menggunakan batas error 0.00001, diperoleh:

exact_value=0.9996645
tol <-  0.00001
err <- 1
m = 1000

while(err>tol){
  set.seed(123)
  hasilCDF<-mc_integral(g,a=0,b=4,m = m)
  err <- abs(hasilCDF-exact_value)
  cat("m=",m,", result=",hasilCDF,", error=",err,"\n",sep = "")
  m=m+1
  if(m==15000){
  break
  }
}
## m=1000, result=1.001826, error=0.002161735
## m=1001, result=1.001721, error=0.00205626
## m=1002, result=1.00079, error=0.00112555
## m=1003, result=1.002007, error=0.002342122
## m=1004, result=1.001017, error=0.001352742
## m=1005, result=1.00003, error=0.0003657328
## m=1006, result=0.99921, error=0.0004544957
## m=1007, result=0.998234, error=0.001430468
## m=1008, result=0.9979907, error=0.001673807
## m=1009, result=1.001692, error=0.002027218
## m=1010, result=1.000933, error=0.001268897
## m=1011, result=1.00014, error=0.0004751195
## m=1012, result=0.9996682, error=3.710246e-06

Jadi, dengan menggunakan Metode Monte Carlo diperoleh nilai \(CDF=0.9996682\) dengan nilai \(error=3.710246e-06\) dan \(m=1012\).

Perbandingan 4 Metode

Berikut merupakan perbandingan hasil perhitungan menggunakan 4 metode di atas dengan pendekatan nilai eksak \(CDF=0.9996645\) dan batas \(error=0.00001\):

Metode <- c("Trapezoidal",
            "Simpson",
            "Gauss Quadrature",
            "Monte Carlo")

Nilai.Eksak<-c(0.9996645,
               0.9996645,
               0.9996645,
               0.9996645)

Nilai_CDF<-c(0.9996745,
             0.9996734,
             0.999662,
             0.9996682)

Nilai_Error<-c(9.98753e-06,
               8.881142e-06,
               2.471959e-06,
               3.710246e-06)

n<-c(732,40,6,"-")

m<-c("-", "-", "-", 1012)

Tabel_perbandingan2 <- data.frame(Metode, Nilai.Eksak, Nilai_CDF, Nilai_Error, n, m)
datatable(Tabel_perbandingan2)

Berdasarkan nilai error dan CDF dari masing-masing metode, dapat dilihat bahwa metode yang menghasilkan nilai paling mendekati nilai eksak dari CDF adalah Metode Gauss Quadrature, karena metode ini menghasilkan nilai error paling kecil dibandingkan dengan metode lainnya. Jadi, dapat disimpulkan bahwa metode integral numerik terbaik untuk menghitung nilai CDF adalah Metode Gauss Quadrature.