################# README file for rcwish ####################
#
# 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.
#
# Copyright January 2000, Phil Everson
# Dept. of Mathematics and Statistics, Swarthmore College. 
#
# 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.
#
#
######################################################################
INSTALLATION:
#############
# There are two pieces to this program, contained in the text files:
# rcwishf.txt, rcwish.txt. 
#
# The first file, rcwishf.txt, contains Fortran code defining subroutines
# called by the SPlus functions in cwish.txt. Save rcwishf.txt to a file
# called rcwish.f, designating it as a Fortran source file. To compile the
# subroutines, type
# 
# f77 -c rcwish.f  (or g77 -c rcwish.f)
#
# at the unix prompt. If you are in a different environment (or have 
# a different version of Fortran), issue the commands necessary to 
# compile subroutines alone. When you finish, you should have created 
# a Fortran object file called rcwish.o.
#
# The second file, cwish.txt, contains SPlus code for the user interface.
# One customization is needed. The call that dynamically loads the Fortran
# subroutine appears in the Splus function rcwish.f as follows:
#
# !!!! 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 !!!!
#
# Edit the path "/home/everson/CWish/" to point to the directory in
# which you have stored rcwish.o. 
# Some versions of SPlus use "dyn.load2" instead of "dyn.load", so you 
# may have to change that as well. If you are running SPlus from other 
# than a unix platform, you should check your manuals about dynamic 
# loading ("loading Fortran and C"). Once you've made these changes, 
# you can source the code into SPlus by typing:

source("cwish.txt")

# from the SPlus prompt, if you happen to store rcwish.txt in your home
# directory, or by specifying the correct path in the source command. for
# example:

source("/home/everson/CWish/rcwish.txt") 

# Once you've done that, you should be good to go. 
######################################################################
#
#
######################################################################
# FUNCTION CALL:
################
# out_rcwish(N, p, df, Sigma=NA, Q=NA, seed=NA, output=0)
#
# REQUIRED INPUTS:
##################
# N:       Number of CWish(df, Sigma; Q) matrices desired.
# p:       Dimension of the desired matrices.
# df:      Degrees of freedom for CWish distribution.
#
# OPTIONAL INPUTS:
##################
# 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.
#
# OUTPUTS:
##########
#  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.
######################################################################
#
#
######################################################################
# EXAMPLE:
##########
#
# Suppose data are observed from the following multivariate Normal
# hierarchical model, with p=5 and k=20:
#
# Level-1:
#  Yi | thetai ~ Np(thetai, I), i = 1,...,k 
#
# Level-2:
#  thetai | A  ~ Np(0, A)
#
# (The Yi and thetai (i=1,...,k) are vectors of length p=5,
#  and A is a 5x5 covariance matrix).
#
# The sufficient statistic S = sum(Yi%*%Yi') has a Wishart
# sampling distribution (given A):  S ~ Wishp(k, I + A).
#
# Suppose (unknown to the statistician) A = I/2: 
# (1/2 times a 5x5 identity matrix)
#
 p_5
 k_20
 I_diag(5)
 A_I/2
#
# Now lets simulate data S using rcwish:
#
 S_rcwish(N=1, p=5, df=k, Sigma=I+A, seed=10000)[,,1]
#
# Leaving Q missing yields unconstrained Wisharts. Note that output 
# from rcwish is an array, even if N=1 (pxpx1).
S
#
#              [,1]       [,2]        [,3]       [,4]          [,5] 
#[1,]  32.887926290  -4.213102  -0.6143759  -1.633734   0.009606083
#[2,]  -4.213102193  26.479641  -7.0343332  -1.644149  -1.013147342
#[3,]  -0.614375932  -7.034333  35.2783991  -2.121467  -5.279387996
#[4,]  -1.633734319  -1.644149  -2.1214672  17.338721   6.025188119
#[5,]   0.009606083  -1.013147  -5.2793880   6.025188  15.597374228
#
# Now the posterior distribution of B = (I + A)^{-1} is Constrained 
# Wishart:  B | S ~ CWishp(k + alpha, S^{-1}; I).
#
# alpha is the prior parameter specified by the statistician (who does
# not know the values of A or B). alpha = 0 is Jeffreys' prior,
# alpha = p+1 is Uniform on B, and alpha = -(p+1) is Uniform on A.
# Let's use alpha = 0 for simplicity. Then the desired CWish df is k=20,
# the scale matrix is S^{-1}, and the constraint matrix is I.
#
Sinv_solve(S)
#
# Now let's simulate 1000 B's from the posterior distribution:
# (this took about 20 seconds on my computer)
#
outB_rcwish(N=1000,p=5,df=k,Sigma=Sinv, Q=I,seed=99999)
#
# Look at the first simulated B:
outB[,,1]
#             [,1]        [,2]         [,3]         [,4]         [,5] 
#[1,]  0.481911436 -0.01181138 -0.130042506 -0.005939071  0.071562764
#[2,] -0.011811381  0.47956631  0.239004185 -0.202304955  0.012466720
#[3,] -0.130042506  0.23900418  0.458766271  0.079449943  0.001351078
#[4,] -0.005939071 -0.20230495  0.079449943  0.779096419 -0.082891542
#[5,]  0.071562764  0.01246672  0.001351078 -0.082891542  0.891353933
#
# The standard output is a (pxpxN) array of the N CWish matrices. To
# save efficiency statistics, set output=1 in the call (default is 0)
#
outB_rcwish(N=1000,p=5,df=k,Sigma=Sinv, Q=I,seed=99999,output=1)
#
# With output=1, the result is a list of length 3. The first component
# is the pxpxN array of CWish matrices. 
outB[[1]][,,1]
#             [,1]        [,2]         [,3]         [,4]         [,5] 
#[1,]  0.481911436 -0.01181138 -0.130042506 -0.005939071  0.071562764
#[2,] -0.011811381  0.47956631  0.239004185 -0.202304955  0.012466720
#[3,] -0.130042506  0.23900418  0.458766271  0.079449943  0.001351078
#[4,] -0.005939071 -0.20230495  0.079449943  0.779096419 -0.082891542
#[5,]  0.071562764  0.01246672  0.001351078 -0.082891542  0.891353933
#
# The second is a 3-vector with the number of matrices begun in order 
# to obtain N acceptable CWish matrices, along the with the number of 
# diagonal and off-diagonal random variables generated.
#
outB[[2]]
# matrices diag rv's off-diag rv's 
#    15959     45389         50245
#
# 15,959 matrices had to be started to obtain 1000 acceptable draws,
# so the overall acceptance rate was 
1000/15959
# [1] 0.06266057
#
# The program uses a stepwise algorithm, so not every matrix begun is 
# completed. If they were, there would have been 
15959*p 
# 79795 diagonal random variables (vs. 45,389 actual) and 
15959*p*(p-1)/2
# 159590 off-diagonal random variables (vs. 50,245 actual). 
#
# The third component of the output is a p-vector recording the 
# number of rejections that occurred at each step of assembly.
outB[[3]]
#  step 1 step 2 step 3 step 4 step 5 
#       0   7971   4366   1761    861
#
# No matrices are rejected at step 1, but 7971 rejections occured
# at step 2 (out of 15,959-1000 = 14,959 total rejections), 4366 at
# step3, and only 861 were completed before rejection.
#
######################################################################