# R/rmf.matrix.gibbs.R In rstiefel: Random Orthonormal Matrix Generation and Optimization on the Stiefel Manifold

#### Documented in rmf.matrix.gibbs

```#' Gibbs Sampling for the Matrix-variate von Mises-Fisher Distribution
#'
#' Simulate a random orthonormal matrix from the matrix von Mises-Fisher
#' distribution using Gibbs sampling.
#'
#'
#' @param M a matrix.
#' @param X the current value of the random orthonormal matrix.
#' @param rscol the number of columns to update simultaneously.
#' @return a new value of the matrix \code{X} obtained by Gibbs sampling.
#' @note This provides one Gibbs scan. The function should be used iteratively.
#' @author Peter Hoff
#' @references Hoff(2009)
#' @examples
#'
#' Z<-matrix(rnorm(10*5),10,5)
#'
#' U<-rmf.matrix(Z)
#' U<-rmf.matrix.gibbs(Z,U)
#'
#'
#' ## The function is currently defined as
#' function (M, X, rscol = NULL)
#' {
#'     if (is.null(rscol)) {
#'         rscol <- max(2, min(round(log(dim(M))), dim(M)))
#'     }
#'     sM <- svd(M)
#'     H <- sM\$u %*% diag(sM\$d)
#'     Y <- X %*% sM\$v
#'     m <- dim(H)
#'     R <- dim(H)
#'     for (iter in 1:round(R/rscol)) {
#'         r <- sample(seq(1, R, length = R), rscol)
#'         N <- NullC(Y[, -r])
#'         y <- rmf.matrix(t(N) %*% H[, r])
#'         Y[, r] <- N %*% y
#'     }
#'     Y %*% t(sM\$v)
#'   }
#'
#' @export rmf.matrix.gibbs
rmf.matrix.gibbs <-
function(M,X,rscol=NULL)
{
#simulate from the matrix mf distribution as described in Hoff(2009)
#this is one Gibbs step, and must be used iteratively

#The number of columns to replace should be small enough so that
#rmf.matrix will be quick, but not so small that you are not taking
#advantage of rmf.matrix. In particular, for square matrices you must
#be replacing two or more columns at a time, or else the chain is
#not irreducible.

if(is.null(rscol)){rscol<-max(2,min( round(log(dim(M))),dim(M))) }

sM<-svd(M)
H<-sM\$u%*%diag(sM\$d)
Y<-X%*%sM\$v

m<-dim(H) ; R<-dim(H)
for(iter in 1:round(R/rscol))
{
r<-sample(seq(1,R,length=R),rscol)
N<-NullC(Y[,-r])
y<-rmf.matrix(t(N)%*%H[,r])
Y[,r]<- N%*%y
}
Y%*%t(sM\$v)
}
```

## Try the rstiefel package in your browser

Any scripts or data that you put into this service are public.

rstiefel documentation built on Feb. 21, 2019, 4:31 p.m.