/* Generates random matrices with a Wishart distribution */
/* Uses the algorithm on pp. 231-232 of W J Kennedy, Jr, */
/* and J E Gentle, Statistical Computing, Dekker 1980. */
/* On input, p is the dimension of the matrix, n the sample size */
/* On output, w is the Wishart matrix and invw its inverse */
/* */
/* Before including this file there must be a statement */
/* #define P defines P as the maximum value of the dimension p */
/* of the symmetric p by p variance-covaraince matrix to be used.*/
/* Often only one value of p is to be used in which case */
/* #define P ... can be followed by int p=P; */
/* */
/* It is necessary to include nrutil.c, jacobi.c and random.c */
/* with this file */
#include "random.h"
#include "nrutil.h"
void wish(int p,int n,float var[P][P],float w[P][P],
float invw[P][P],long *idum)
{
int i,j,k,nrot;
float gamdev(int ia, long *idum);
float gasdev(long *idum);
float *d,*y,**v,**evec,**a,**z,**b,**s,
**invs,**intermed,temp;
d=vector(1,P);
y=vector(1,P);
evec=matrix(1,P,1,P);
a=matrix(1,P,1,P);
z=matrix(1,P,1,P);
b=matrix(1,P,1,P);
s=matrix(1,P,1,P);
invs=matrix(1,P,1,P);
intermed=matrix(1,P,1,P);
v=convert_matrix(&var[0][0],1,p,1,p);
/* Ensure symmetry */
for (i=1;i <= p;i++)
for (j=1;j <= p;j++)
if (i < j) v[i][j]=v[j][i];
/* Step 1. Determine a such that v = a*+Transpose(a) */
jacobi(v,p,d,evec,&nrot);
/* d consists of eigenvalues, evec of eigenvectors */
for (i=1;i <= p;i++)
for (j=1;j <= p;j++)
a[i][j]=evec[i][j]*sqrt(d[j]);
/* Step 2. Generate z[i][j] ~ N(0,1) */
/* for i <= j = 1, 2, ..., p */
for (j=1;j <= p;j++)
for (i=1;i <= j;i++)
z[i][j]=gasdev(idum);
/* Step 3. Generate y[i] ~ chi^2_{n-i} */
/* for i = 1, 2, ..., p */
for (i=1;i <= p;i++)
y[i]=2*gamdev((n-i)/2,idum);
/* Step 4. Find b as in W J Kennedy, Jr, */
/* & J E Gentle, p. 231 */
b[1][1]=y[1];
for (j=2;j <= p;j++) {
b[j][j]=y[j];
for (i=1;i < j;i++)
b[j][j]+=z[i][j]*z[i][j];
b[1][j]=z[1][j]*sqrt(y[i]);
for (i=1;i < j;i++) {
b[i][j]=z[i][j]*sqrt(y[i]);
for (k=1;k < i;k++)
b[i][j]+=z[k][i]*z[k][j];
}
for (i=1;i <= j;i++)
b[i][j]=b[j][i];
}
/* Step 5. Find s = a*+b*+Transpose(a) */
/* First intermed = b*+Transpose(a) */
for (i=1;i <= p;i++)
for (j=1;j <= p;j++) {
intermed[i][j]=0;
for (k=1;k <= p;k++)
intermed[i][j]+=b[i][k]*a[j][k];
}
/* Now s = a*+intermed */
for (i=1;i <= p;i++)
for (j=1;j <= p;j++) {
s[i][j]=0;
for (k=1;k <= p;k++)
s[i][j]+=a[i][k]*intermed[k][j];
}
/* Now ensure symmetry */
for (i=1;i <= p;i++)
for (j=1;j < i;j++) {
temp=(s[i][j]+s[j][i])/2;
s[i][j]=temp;
s[j][i]=temp;
}
/* Now s is a Wishart matrix (aside */
/* from the factor 1/n) with n-p d.f. */
/* (W J Kennedy, Jr, & J E Gentle, p. 232) */
/* Step 6. Find inverse invs of s */
/* First determine a such that */
/* s = a*+Transpose(a) */
jacobi(s,p,d,evec,&nrot);
/* Resymmetrize s */
for (i=1;i <= p;i++)
for (j=1;j < i;j++)
s[j][i]=s[i][j];
/* d consists of eigenvalues, evec of */
/* eigenvectors. From this, find inverse */
for (i=1;i <= p;i++)
for (j=1;j <= p;j++) {
invs[i][j]=0;
for (k=1;k <= p;k++)
invs[i][j]+=evec[i][k]*(1/d[k])*evec[j][k];
}
for (i=0;i < p;i++)
for (j=0;j < p;j++)
w[i][j]=s[i+1][j+1];
for (i=0;i < p;i++)
for (j=0;j < p;j++)
invw[i][j]=invs[i+1][j+1];
free_convert_matrix(v,1,p,1,p);
free_matrix(evec,1,P,1,P);
free_matrix(a,1,P,1,P);
free_matrix(z,1,P,1,P);
free_matrix(b,1,P,1,P);
free_matrix(s,1,P,1,P);
free_matrix(invs,1,P,1,P);
free_matrix(intermed,1,P,1,P);
free_vector(d,1,P);
free_vector(y,1,P);
}
/* DISCLAIMER OF WARRANTY */
/* THE PROGRAMS ACCESSED BY THIS ROUTINE (AND ON THE ORIGINAL */
/* DISKETTE) ARE PROVIDED "AS IS" WITHOUT WARRANTY OF ANY KIND. */
/* WE MAKE NO WARRANTIES, EXPRESS OR IMPLIED, THAT THEY ARE FREE OF */
/* ERROR, OR ARE CONSISTENT WITH ANY PARTICULAR STANDARD OF */
/* MERCHANTABILITY, OR THAT THEY WILL MEET YOUR REQUIREMENTS FOR ANY */
/* PARTICULAR APPLICATION. THEY SHOULD NOT BE RELIED ON FOR SOLVING */
/* A PROBLEM WHOSE INCORRECT SOLUTION COULD RESULT IN INJURY TO A */
/* PERSON OR LOSS OF PROPERTY. IF YOU DO USE THEM IN SUCH A MANNER, */
/* IT IS AT YOUR OWN RISK. THE AUTHORS AND PUBLISHER DISCLAIM ALL */
/* LIABILITY FOR DIRECT, INDIRECT, OR CONSEQUENTIAL DAMAGES */
/* RESULTING FROM YOUR USE OF THE PROGRAMS. */
/* */