function [theta1, theta2]=gau(k,amu,h);
   global a1;
   a1=(chol(inv(h)))';  % a1*a1'=h
   v1=a1*randn(k,1);
   theta1=amu+v1';
   theta2=amu-v1';   
end;