Num Fun!

Numerical Analysis in R - Sorting, Searching, Root Finding and Integration

CSJP

Sorting Begins …

Selection Sort

set.seed(1234)
a = sample(10)
a
##  [1]  2  6  5  8  9  4  1  7 10  3
selsort = function(x) {
    for (i in 1:(length(x) - 1)) {
        min = i
        for (j in (i + 1):length(x)) if (x[j] < x[min]) 
            min = j  # find the ith smallest number 
        t = x[min]  # make exchanges btw ith smallest one and ith position element 
        x[min] = x[i]
        x[i] = t
    }
    x
}
selsort(a)
##  [1]  1  2  3  4  5  6  7  8  9 10

Searching Steps in …

set.seed(1234)
v = sample(1:26, 1)
v
## [1] 3
binsearch = function(x) {
    p = 1
    r = length(LETTERS[1:26])
    if (x <= 0 | x > r) 
        return(print("No Such LETTER!"))
    while (r >= p) {
        h = floor((p + r)/2)  # divide and conquer
        if (x < h) 
            r = h - 1  # fall in lower half
 else p = h + 1  # fall in upper half
        if (x == h) 
            return(LETTERS[h])  # find it
    }
}

binsearch(v)
## [1] "C"

v = 30
binsearch(v)
## [1] "No Such LETTER!"

v = 26
binsearch(v)
## [1] "Z"

v = 14
binsearch(v)
## [1] "N"

Finding Roots …

Plot Function f(x)

# Plot Function f(x)
f = function(x) x^3 - 9 * x
curve(f, -5, 5, ylab = expression(f(x) == x^3 - 9 %*% x), main = expression(paste("Plot of ", 
    f(x) == x^3 - 9 %*% x)))
abline(h = 0, v = 0)
grid()
points(c(-3, 0, 3), c(0, 0, 0), pch = 20)
text(c(-3, -0.35, 3), c(8, -8, -8), c("(-3, 0)", "(0, 0)", "(3, 0)"))

Root Finding (1) …

Bisect Method

# Bisect method
bisect = function(f, x0, x1, tol, itmax) {
    if (sign(f(x0)) == sign(f(x1))) 
        return(print("Wrong Initial Conditions!"))
    for (i in 1:itmax) {
        x = (x0 + x1)/2
        if (abs(f(x)) < tol) 
            return(x)
        if (sign(f(x)) == sign(f(x0))) 
            x0 = x else x1 = x
    }
    return(print("Root Not Found!"))
}
bisect(f, -4, -2, 1e-05, 100)
## [1] -3
bisect(f, 2, 4, 1e-05, 100)
## [1] 3

Root Finding (2) …

Secant Method

# Secant method
f = function(x) x^3 - 9 * x
secant = function(f, x0, x1, tol, itmax) {
    for (i in 1:itmax) {
        x = x1 - f(x1) * (x1 - x0)/(f(x1) - f(x0))
        if (abs(f(x)) < tol) 
            return(x)
        x0 = x1
        x1 = x
    }
    return(print("Root Not Found!"))
}
secant(f, -4, -2, 1e-05, 100)
## [1] -3
secant(f, 2, 4, 1e-05, 100)
## [1] 3

Root Finding (3) …

Newton Method

# Newton method
f = expression(x^3 - 9 * x)
fprime = D(f, "x")
newton = function(f, fprime, tol, itmax) {
    for (i in 1:itmax) {
        x = x - eval(f)/eval(fprime)
        if (abs(eval(f)) < tol) 
            return(x)
    }
    return(print("Root Not Found!"))
}
x = -2
newton(f, fprime, 1e-05, 100)
## [1] -3
x = 2
newton(f, fprime, 1e-05, 100)
## [1] 3

Finally, Integration Comes …

Using R function

# Using R function
f = function(x) 1/x
fint = integrate(f, 1, 2)
eval(fint)
## 0.6931 with absolute error < 7.7e-15

Simple Quadrature Method

# simple quadrature method
f = function(x) 1/x
intrect = function(f, a, b, N) {
    r = 0
    w = (b - a)/N
    for (i in 1:N) r = r + w * f(a - w/2 + i * w)
    return(r)
}
intrect(f, 1, 2, 10)
## [1] 0.6928
intrect(f, 1, 2, 1000)
## [1] 0.6931

my sessionInfo …

print(sessionInfo(), locale = FALSE)
## R version 3.0.2 (2013-09-25)
## Platform: i386-w64-mingw32/i386 (32-bit)
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] cluster_1.14.4     foreach_1.4.1      plyr_1.8          
##  [4] reshape2_1.2.2     caret_6.0-21       lattice_0.20-24   
##  [7] randomForest_4.6-7 e1071_1.6-2        class_7.3-9       
## [10] nnet_7.3-7         ada_2.0-3          rpart_4.1-4       
## [13] gtools_3.2.1       gdata_2.13.2       caTools_1.16      
## [16] KernSmooth_2.23-10 MASS_7.3-29        ROCR_1.0-5        
## [19] gplots_2.12.1      ggplot2_0.9.3.1    scales_0.2.3      
## [22] RColorBrewer_1.0-5 colorspace_1.2-4   vcd_1.3-1         
## [25] boot_1.3-9         knitr_1.5         
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-6    car_2.0-19      codetools_0.2-8 dichromat_2.0-0
##  [5] digest_0.6.4    evaluate_0.5.1  formatR_0.10    gtable_0.1.2   
##  [9] iterators_1.0.6 labeling_0.2    munsell_0.4.2   proto_0.3-10   
## [13] stringr_0.6.2   tools_3.0.2