## Wednesday, August 20, 2008

### More stats sites

Jim Albert has a blog based on his book Bayesian Computation in R. Unfortunately, it seems they haven't been updating it too much recently but it looks like what's there is pretty good.

## Monday, August 18, 2008

### Matrix square roots in R

There are a couple of ways to do matrix square roots in R. The first is through eigendecomposition:

# a is the matrix for which we want to square root:
a.eig <- eigen(a)
a.sqrt <- a.eig\$vectors %*% diag(sqrt(a.eig\$values)) %*% solve(a.eig\$vectors)

Of course, this only works for positive semidefinite matrices that are diagonalizable (no normal Jordan forms, etc.)

Another way is to use the Denman-Beavers algorithm:

denman.beavers <- function(mat,maxit=50) {
stopifnot(nrow(mat) == ncol(mat))
niter <- 0
y <- mat
z <- diag(rep(1,nrow(mat)))
for (niter in 1:maxit) {
y.temp <- 0.5*(y+solve(z))
z <- 0.5*(z+solve(y))
y <- y.temp
}
return(list(sqrt=y,sqrt.inv=z))
}

This can probably be improved with better convergence criteria, but so far it works well and gives both the square root and inverse square root, both of which are required for many applications.

## Saturday, August 16, 2008

### Neat, someone else in the niche

On biostatistics and clinical trial. I'm surprised I didn't discover him earlier.