## 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.

# 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.

## Wednesday, August 13, 2008

## 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.

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

Subscribe to:
Posts (Atom)