function p = invgauss(x, mu, D)
    if abs(x) < 1e-8
        p = zeros(size(x))
    else
        if x < 0
            x = -x;
        end
        p = exp(-(x-mu).^2./4.0./D./x./mu./mu)./sqrt(4.*pi.*D.*x.^3.0);
    end
    p(p<1e-16) = 1e-16;
end