# rbing.matrix.gibbs: Gibbs Sampling for the Matrix-variate Bingham Distribution In rstiefel: Random orthonormal matrix generation on the Stiefel manifold

## Description

Simulate a random orthonormal matrix from the Bingham distribution using Gibbs sampling.

## Usage

 `1` ```rbing.matrix.gibbs(A, B, X) ```

## Arguments

 `A` a symmetric matrix. `B` a diagonal matrix with decreasing entries. `X` the current value of the random orthonormal matrix.

## Value

a new value of the matrix `X` obtained by Gibbs sampling.

## Note

This provides one Gibbs scan. The function should be used iteratively.

Peter Hoff

Hoff(2009)

## Examples

 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28``` ```Z<-matrix(rnorm(10*5),10,5) ; A<-t(Z)%*%Z B<-diag(sort(rexp(5),decreasing=TRUE)) U<-rbing.Op(A,B) U<-rbing.matrix.gibbs(A,B,U) ## The function is currently defined as function (A, B, X) { m <- dim(X)[1] R <- dim(X)[2] if (m > R) { for (r in sample(seq(1, R, length = R))) { N <- NullC(X[, -r]) An <- B[r, r] * t(N) %*% (A) %*% N X[, r] <- N %*% rbing.vector.gibbs(An, t(N) %*% X[, r]) } } if (m == R) { for (s in seq(1, R, length = R)) { r <- sort(sample(seq(1, R, length = R), 2)) N <- NullC(X[, -r]) An <- t(N) %*% A %*% N X[, r] <- N %*% rbing.Op(An, B[r, r]) } } X } ```

### Example output

```function (A, B, X)
{
m <- dim(X)[1]
R <- dim(X)[2]
if (m > R) {
for (r in sample(seq(1, R, length = R))) {
N <- NullC(X[, -r])
An <- B[r, r] * t(N) %*% (A) %*% N
X[, r] <- N %*% rbing.vector.gibbs(An, t(N) %*% X[,
r])
}
}
if (m == R) {
for (s in seq(1, R, length = R)) {
r <- sort(sample(seq(1, R, length = R), 2))
N <- NullC(X[, -r])
An <- t(N) %*% A %*% N
X[, r] <- N %*% rbing.Op(An, B[r, r])
}
}
X
}
```

rstiefel documentation built on May 29, 2017, 11:58 a.m.