function x=rand_nrig(n,g)
% count 'n' natural, shape 'g' positive real
% generate standard normal reciprocal inverse gaussian variates
% by mixing normals whose variances are mean-one reciprocal
% inverse gaussian variates generated by the method of
% Michael-Schucany-Haas[1976]
% JDodson 20151009
    y=randn(n,1);
    p=(1+1/g)/(2*g)*(2*g+y.^2-abs(y).*sqrt(4*g+y.^2));
    ind=find(rand(n,1)>1 ./(1+p/(1+1/g)));
    p(ind)=(1+1/g)^2 ./p(ind);
    x=randn(n,1)./sqrt(p);
end