gauk[k_, amu_, h_, x_] := - (0.5)*(x-amu).h.(x-amu);
gaun[k_, h_] := -(0.5)*k*Log[2*N[Pi]] + (0.5)*Log[Abs[Det[h]]]//N;
gau[k_, amu_, h_] := (
ndist=NormalDistribution[0,1];
x=Table[Random[ndist],{k}];
var=Inverse[h];
a1=Eigenvectors[N[var]].DiagonalMatrix[(Eigenvalues[N[var]])^(1/2)];
(* a1.a1'=var *)
v1=a1.x;
theta1=amu+v1;
theta2=amu-v1;
)
gaua[k_, amu_, h_] := (
x=Table[Random[ndist],{k}];
v1=a1.x;
theta1=amu+v1;
theta2=amu-v1;
)