# rbing.Op: Simulate a 'p*p' Orthogonal Random Matrix In rstiefel: Random Orthonormal Matrix Generation and Optimization on the Stiefel Manifold

## Description

Simulate a `p*p` random orthogonal matrix from the Bingham distribution using a rejection sampler.

## Usage

 `1` ```rbing.Op(A, B) ```

## Arguments

 `A` a symmetric matrix. `B` a diagonal matrix with decreasing entries.

## Value

A random pxp orthogonal matrix simulated from the Bingham distribution.

## Note

This only works for small matrices, otherwise the sampler will reject too frequently to be useful.

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 29 30 31 32 33 34 35 36 37 38 39 40 41``` ```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) { b <- diag(B) bmx <- max(b) bmn <- min(b) if(bmx>bmn) { A <- A * (bmx - bmn) b <- (b - bmn)/(bmx - bmn) vlA <- eigen(A)\$val diag(A) <- diag(A) - vlA[1] vlA <- eigen(A)\$val nu <- max(dim(A)[1] + 1, round(-vlA[length(vlA)])) del <- nu/2 M <- solve(diag(del, nrow = dim(A)[1]) - A)/2 rej <- TRUE cholM <- chol(M) nrej <- 0 while (rej) { Z <- matrix(rnorm(nu * dim(M)[1]), nrow = nu, ncol = dim(M)[1]) Y <- Z %*% cholM tmp <- eigen(t(Y) %*% Y) U <- tmp\$vec %*% diag((-1)^rbinom(dim(A)[1], 1, 0.5)) L <- diag(tmp\$val) D <- diag(b) - L lrr <- sum(diag((D %*% t(U) %*% A %*% U))) - sum(-sort(diag(-D)) * vlA) rej <- (log(runif(1)) > lrr) nrej <- nrej + 1 } } if(bmx==bmn) { U<-rustiefel(dim(A)[1],dim(A)[1]) } U } ```

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