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

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.

Wednesday, August 13, 2008


Does a statistician need to twitter? Who knows. But why not try.

Monday, August 4, 2008

Lan-DeMets in R is easier

I found this little gem today. Makes simulating group sequential designs a lot easier, since I don't have to run ld98 every time or something silly like that.

In other news, I got some cool toys from Google at JSM.