Numerical Analysis talk at Statistical Programming DC

I gave my talk tonight on numerical analysis in R and the talk is already up on YouTube.

There’s an interesting set of questions at the end about resources for introducing new users to R. The slides are not terribly visible, so you can see them here:

And this is the R script used for demonstrations:

## Example of vector arithmetic
u <- c(1, 2, 3); v <- c(8, 4, 2); x <- 7
u + x
u + v

## Element recycling
u + c(1, 9)

## Matrix arithmetic
A <- matrix(1:9, 3)
A + 1
A + u
A + c(1, 2)

## Diagonalization
(B <- matrix(1:6, 3))
diag(A)
diag(B)

## Matrix multiplication
A %*% B

## Cross product 
crossprod(A, B)

## Outer product
outer(u, v)

## dot product
u %*% u 

## Matrix solution / RREF
(A <- matrix(c(2, 3, 1, 1, 2, -5, -1, -2, 4), 3))
(b <- c(1, 1, 3))
solve(A, b)

## LU Decomposition
library(Matrix)
lu(A)

## Linear interpolation
x <- 1:10
y <- log(x)
plot(x, y, main = "approxfun demo")
par(new = TRUE)
flog <- approxfun(x, y)
plot(flog, 1, 10, col = "blue")

## Spline
x <- seq(-10, 10, 4)
y <- x^2 + 5 * x - 3
plot(x, y, main = "splinefun demo")
par(new = TRUE)
linfun <- approxfun(x, y)
plot(linfun, -10, 10, col = "blue")
par(new = TRUE)
linfun <- splinefun(x, y)
plot(linfun, -10, 10, col = "green")

## Polynomial interpolation
polyinterp <- function(x, y) {
  if(length(x) != length(y))
    stop("Length of x and y vectors must be the same")
  
  n <- length(x) - 1
  vandermonde <- rep(1, length(x))
  for(i in 1:n) {
    xi <- x^i
    vandermonde <- cbind(xi, vandermonde)
  }
  beta <- solve(vandermonde, y)
  
  names(beta) <- NULL
  return(rev(beta))
}

x <- c(1, 2, 3)
y <- x^2 + 5 * x - 3
z <- polyinterp(x, y)

rhorner(1, z)

## Integration
integrate(log, 2, 6)
integrate(flog, 2, 6)

## Root Finding
## x^2 + x - 6
f <- function(x) { (x + 3) * (x - 2) }             
plot(f, -5, 5, col = "blue")
abline(h = 0, v = 0, col = "black", lty = "dotted") 
uniroot(f, c(-5, 0))

z <- c(-6, 1, 1)
polyroot(z)

## x^2 + x + 6
z <- c(6, 1, 1)
polyroot(z)

## Horner's method
rhorner <- function(x, betas) {
  n <- length(betas)
  
  if(n == 1)
    return(betas)
  
  return(betas[1] + x * rhorner(x, betas[2:n]))
}

rhorner(z, 5)

Many thanks to Marck Vaisman for putting this together and Casey Patrick Driscoll for putting it on YouTube!