proc gau(k,amu,h);
   a1=(chol(invpd(h)))';  /* a1*a1'=h */
   v1=a1*rndn(k,1); 
   theta1=amu+v1';
   theta2=amu-v1';
   retp(" ");
endp;