1

I am trying to draw random samples from some distribution as follows:enter image description here

my code runs but the numbers look strange. so I am not sure what went wrong, maybe some operators. The elements are extremely large. my attempt:

C_hat=(((x`)*x)**(-1))*((x`)*z);
S=((z-x*c_hat)`)*((z-x*c_hat));

*draw sigma;
sigma = shape(RandWishart(1, 513 - 3 - 2,s**(-1)),4,4);
*draw vec(c);
vec_c_hat= colvec(c_hat`); *vectorization of c_hat;
call randseed(4321);
vec_c = RandNormal(1,vec_c_hat,(sigma`)@(((x`)*x)**(-1)));
c = shape(vec_c,4,4);
print c;
duckman
  • 687
  • 1
  • 15
  • 30

1 Answers1

0

Since you haven't provided data or a reference, it is difficult to guess whether your "strange" and "extremely large" numbers are correct. However, the program looks mostly correct, so check your data.

A minor problem with your program is that you are using the SHAPE function to reshape the vec_c vector into a matrix. You should be using the SHAPECOL function (or transpose the result).

The following program uses the Sashelp.Cars data, which is distributed with SAS, to initialize the X and Z matrices. The program computes a random matrix C which is close to the inverse crossproduct matrix for the data. I've also added some intermediate computations and comments. This version works as expected on the Sashelp.Cars data:

proc iml;
use sashelp.cars;
read all var {weight wheelbase enginesize horsepower} into X;
read all var {mpg_city mpg_highway} into Z;
close;

*normal equations and covariance;
xpx = x`*x;
invxpx = inv(xpx);
C_hat = invxpx*(x`*z);
r = z-x*c_hat;
S = r`*r;

*draw sigma;
call randseed(4321);
DF = nrow(X)-ncol(X)-2;
W = RandWishart(1, DF, inv(S));  /* 1 x (p*p) row vector */
sigma = shape(W, sqrt(ncol(W))); /* reshape to p x p matrix */
*draw vec(c);
vec_c_hat = colvec(c_hat`);      /* stack columns of c_hat */
vec_c = RandNormal(1, vec_c_hat, sigma@invxpx);
c = shapecol(vec_c, nrow(C_hat), ncol(C_hat));   /* reshape w/ SHAPECOL */
print C_hat, c;
Rick
  • 1,210
  • 6
  • 11