### S-PLUS CONSTRAINED WISHART RANDOM MATRIX GENERATOR ### # Copyright January 2000, Phil Everson # Dept. of Mathematics and Statistics, Swarthmore College. # rcwish_function(N, p, df, Sigma=NA, Q=NA, seed=NA, output=0){ # # rcwish is an Splus function which calls Fortran subroutines # (contained in the Fortran object file rcwish.o) to generate # random draws from the Constrained Wishart distribution. # See Everson & Morris (2000), JCGS for details. # if(p>1) # # !!!! IMPORTANT !!!! IMPORTANT !!!! IMPORTANT !!!! IMPORTANT !!!! # dyn.load("/home/everson/CWish/rcwish.o") # ---> ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ <--- # This PATH MUST BE CHANGED to point to where rcwish.o is loacated. # Fortran routines are unneccessary for univariate generation. # # !!!! IMPORTANT !!!! IMPORTANT !!!! IMPORTANT !!!! IMPORTANT !!!! # # This program is free and may be redistributed as long as # the copyright is preserved. # # This program should be cited if used in published research. # # The author will appreciate notification of any comments or # errors by emailing peverso1@swarthmore.edu # #Reference: # Everson & Morris (2000), "Simulation from Wishart distributions # with eigenvalue constraints", Journal of Computational and # Graphical Statistics (to appear). # # Fortran subroutines for Uniform random number generation and for # evaluating the gamma and incomplete gamma functions were taken # from: # # Press, W. H., Teukolsky, S. A., Vetterling, W. T., # Flannery, B. P. (1992), "Numerical Recipes in Fortran" # (2nd ed.), Cambridge University Press. # # # Inputs: # N: Number of CWish(df, Sigma; Q) matrices desired. # p: Dimension of the desired matrices. # df: Degrees of freedom for CWish distribution. # Sigma: Symmetric positive-definite pxp scale matrix. If left # missing, Sigma is set to I, the pxp identity matrix. # Q: Symmetric positive-definite pxp constraint matrix. # Q may also be a p-vector representing a diagonal matrix, # a scalar representing a matrix proportional to I, or left # missing, which yields unconstrained Wisharts. # seed: Integer seed for random number generator. # # output: Controls program output. # # output = 0: Returns a pxpxN array Wa, where Wa[,,i] gives the # ith of N pxp CWish matrices. # output = 1: Returns a list of length 3: 1) Wa; 2) nvec - a # 3-vector containing the number of matrices attempted # and the number of diagonal and off-diagonal random # variables needed; 3) nrej - a p-vector containing the # number of attempts that were rejected at each step. # if(missing(seed)) seed_rpois(1,1e5*rgamma(1,1e4)) # if(p==1){ # # UNIVARIATE PROBLEM # set.seed(seed) nvec_c(N,N,0) nrej_0 if(missing(Sigma)){ Sigma_1 } else{ if(length(Sigma)>1||Sigma<0) stop("Sigma must be a positive scalar when p=1.") } if(missing(Q)){ # # Generate (scaled) unconstrained Chi-squares. # Wa_rchisq(N,df)*Sigma } else{ # # Generate (scaled) Constrained Chi-squares. # if(length(Q) > 1||Q<0) stop("Q must be a positive scalar when p=1.") d_Q/Sigma Wa_rcchisq(N,df,d)*Sigma } } # End Univariate problem (end if(p==1){...}). if(p>1){ # # MULTIVARIATE PROBLEM # standard_0 if(missing(Sigma)){ standard_1 Sigma_diag(p) } else{ eigS_checkSig(Sigma,p) if(max(abs(Sigma-diag(p)))<1e-6) standard_1 } # End Sigma checks. standard=1 indicates Sigma = I. if(missing(Q)){ # # Prepare to generate unconstrained Wisharts. # d_rep(0,p) pd_rep(1,p) if(standard==0){ C_eigS[[2]]%*%diag(sqrt(eigS[[1]]))%*%t(eigS[[2]]) # C%*%C=Sigma } } # else{ # # Prepare to generate Constrained Wisharts. # eigQ_checkQ(Q,p) rtQ_eigQ[[2]]%*%diag(sqrt(eigQ[[1]]))%*%t(eigQ[[2]]) rtQinv_eigQ[[2]]%*%diag(1/sqrt(eigQ[[1]]))%*%t(eigQ[[2]]) # rtQ is a pxp symmetric square root matrix for Q. # rtQinv is a pxp symmetric inverse square root matrix for Q. if(standard==0){ eig_eigen(rtQinv%*%Sigma%*%rtQinv,T) d_1/eig[[1]] Gam_eig[[2]] C_rtQ%*%Gam%*%diag(1/sqrt(d)) } else{ d_eigQ[seq(p,1,-1)] } pd_pchisq(d,df) } # end if(missing(Q)){...}; else{...}). # outU_rscwish(N,p,df,d,pd,seed) nvec_outU[[2]] nrej_outU[[3]] if(standard==1){ Wa_outU[[1]] } else{ Wa_mammult(C,outU[[1]],t(C)) } } # End Multivariate problem (end if(p > 1){...}). # if(output==0) out_Wa if(output==1){ names(nvec)_c("matrices", "diag rv's", "off-diag rv's") names(nrej)_paste("step",1:p) out_list(Wa, nvec, nrej) names(out)_c("CWish", "counts", "rejections") } out } rcchisq_function(N, df, d=NA){ # Returns N random draws from a Chi-Square distribution with df # degrees of freedom, and constrained to be less than d (scalar). if(missing(d)) # unconstrained out_rchisq(N,df) else{ out_qchisq(runif(N)*pchisq(d,df),df) out[out>d]_d } out} rscwish_function(N,p,df,d,pd,seed){ # Calls Fortran subroutine "rscwish" to generate N pxp # standard constrained Wishart matrices with df degrees of # freedom and a p-vector of constraints d. pd is the p-vector # of the Chi-square(df) CDF evaluated at d. If all elements of # d are 0 and all elements of pd are 1, then rscwish returns N # unconstrained Wisharts. # Ua_array(0.0,c(p,p,N)) Tmat_matrix(0.0,p,p) Z_rep(0,p) chitab_c(qchisq(seq(.001,.999,.001), df)) # chitab is a table of 999 Chi-square(df) quantiles # to pass to Fortran. out_.Fortran("rscwish", Ua=Ua, as.integer(N), as.integer(p), as.double(df), as.double(d), as.double(pd), as.double(chitab), as.integer(seed), as.double(Tmat), # Tmat as.double(Tmat), # Deltinv as.double(Tmat), # H11 as.double(Z), # H12 as.double(Z), # Zt as.double(Z), # Zts as.double(Z), # tempp as.integer(0), # mcnt as.integer(0), # drvcnt as.integer(0), # orvcnt as.integer(Z) # rejcnt ) list(out[[1]],c(out[[16]],out[[17]],out[[18]]),c(out[[19]])) } mammult_function(M,A,tM){ # Calls Fortran subroutine "mammult" to pre-multiply by M # (pxp) and post-multiply by tM (pxp) each of the N pxp # elements of A (pxpxN). # pp_dim(A) p_pp[1] N_pp[3] MAM_A*0.0 out_.Fortran("mammult", as.double(M), as.double(tM), as.double(A), as.integer(p), as.integer(N), as.double(0*M), # U as.double(0*M), # V MAM=MAM ) out$MAM } checkSig_function(Sigma,p){ # For multivatiate problems, check that Sigma is symmetric, # and non-negative definite. Returns the p eigenvalues # (sorted largest to smallest) and corresponding eigenvectors. # pp_dim(Sigma) if(length(pp)!=2){ stop("Sigma must be a symmetric, non-neg def, pxp matrix.") } else{ if(pp[1]!=p|pp[2]!=p) stop("Sigma must be a symmetric, non-neg def, pxp matrix.") } if(max(abs(t(Sigma)-Sigma))>.001) stop("Sigma must be a symmetric, non-neg def, pxp matrix.") eig_eigen(Sigma,T) if(min(eig[[1]])<=0) stop("Sigma must be a symmetric, non-neg def, pxp matrix.") eig } checkQ_function(Q,p){ # For constrained multivatiate problems, check the form of Q. # Returns the p eigenvalues (sorted largest to smallest) # and corresponding eigenvectors. if(length(Q)==1){ if(Q<0) stop("Q must be a positive scalar, p-vector, or pxp matrix.") eign_rep(Q,p) rtQ_diag(p)*sqrt(Q) rtQinv_diag(p)/sqrt(Q) } else{ if(length(Q)==p){ eign_sort(Q)[seq(p,1,-1)] if(eign[p]<0){ stop("Q must be a positive scalar, p-vector, or pxp matrix.") } rtQinv_diag(1/sqrt(Q)) } else{ pp_dim(Q) if(length(pp)!=2){ stop("Q must be a positive scalar, p-vector, or pxp matrix.") } else{ if(pp[1]!=p|pp[2]!=p){ stop("Q must be a positive scalar, p-vector, or pxp matrix.") } else{ if(max(abs(t(Q)-Q))>.001) stop("Q must be a symmetric.") eig_eigen(Q,T) if(eig[[1]][p]<0) stop("Q must be non-negative definite.") } } } } eig } rwish_function(N, p, df, Sigma=NA,seed=NA){ # Returns N random p-dimensional Wishart matrix with df degrees # of freedom and scale matrix Sigma (default is identity). ###### This function does *not* call Fortran ######### # if(missing(seed)) seed_.Random.seed[1] # set.seed(seed) if(p==1){ # # UNIVARIATE PROBLEM # if(missing(Sigma)){ Sigma_1 } else{ if(length(Sigma)>1||Sigma<0) stop("Sigma must be a positive scalar when p=1.") } Wa_rchisq(N,df)*Sigma } # End Univariate problem (end if(p==1){...}). if(p>1){ # # MULTIVARIATE PROBLEM # standard_0 if(missing(Sigma)){ standard_1 Sigma_rtSig_diag(p) eign_rep(1,p) } else{ eigS_checkSig(Sigma,p) if(max(abs(Sigma-diag(p)))<1e-6) standard_1 Sigma_(Sigma+t(Sigma))/2 rtSig_eigS[[2]]%*%diag(eigS[[1]]^(0.5))%*%t(eigS[[2]]) } Wa_array(0,c(p,p,N)) for(index in 1:N){ u_rchisq(p, (df-c(1:p)+1) ) Tmat_diag(sqrt(u)) for(i in 2:p){ for(j in 1:(i-1)){ Tmat[i,j]_rnorm(1) } } Tmat_rtSig%*%Tmat Wa[,,index]_Tmat%*%t(Tmat) } } # End Multivariate problem (end if(p>1){...}). Wa }