function zz = rinvnorm(n, mu, lambda)
  zz = zeros(1,n);
  v = randn(1,n);
  y = v.*v;
  %keyboard
  x = mu+mu.*mu.*y./lambda/2-mu/2./lambda.*sqrt(4*mu.*lambda.*y+mu.*mu.*y.*y);
  z = rand(1,n);
  ind = z <= (mu./(mu+x));
  zz(ind) = x(ind);
  tmp = mu.*mu./x;
  zz(~ind) = tmp(~ind);