Biostatistics, clinical trial design, critical thinking about drugs and healthcare, skepticism, the scientific process.
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)