Numerical Analysis talk at Statistical Programming DC | James Howard Numerical Analysis talk at Statistical Programming DC | James Howard

James Howard A Mathematician, a Different Kind of Mathematician, and a Statistician

image representing a theme in this article

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:

Numerical Analysis in R from James Howard

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 = &quot;approxfun demo&quot;)
par(new = TRUE)
flog <- approxfun(x, y)
plot(flog, 1, 10, col = &quot;blue&quot;)

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

## Polynomial interpolation
polyinterp <- function(x, y) {
  if(length(x) != length(y))
    stop(&quot;Length of x and y vectors must be the same&quot;)
  
  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 = &quot;blue&quot;)
abline(h = 0, v = 0, col = &quot;black&quot;, lty = &quot;dotted&quot;) 
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!