using SpecialFunctions, QuadGK, Gadfly

Normal random variable¶

"standard Normal density"
function Normal_dens(x)
    return exp(-x^2/2)/sqrt(2π)
end;

"standard Normal distribution"
function Normal_dist(x)
    return erfc(x/-sqrt(2))/2
end;

plot([Normal_dens,Normal_dist],-2,2)

"standard Normal quantile"
function Normal_quan(p)
    return -sqrt(2)erfcinv(2p)
end;

plot(Normal_quan,.01,.99)

symmetric Normal / reciprocal inverse Gaussian (NRIG) random variable¶

"standard NRIG density"
function NRIG_dens(x,g)
    return exp(g)sqrt(1+g)besselk(0,sqrt(g^2+(1+g)*x^2))/π
end;

"standard NRIG distribution"
function NRIG_dist(x,g)
    return 1/2+quadgk(z->NRIG_dens(z,g),0,x)[1]
end;

plot([x->NRIG_dens(x,1),x->NRIG_dist(x,1)],-2,2)

"standard NRIG quantile"
function NRIG_quan(p,g)
# Newton-Raphson method to invert NRIG_dist
    x0=Normal_quan(p) # initial guess
    x1=NaN
    for n=1:20
        f=NRIG_dens(x0,g)
        x1=(p-1/2+quadgk(x->f-NRIG_dens(x,g),0,x0)[1])/f
        if abs(x1-x0)<1e-12
            break
        end
        x0=x1
    end
    return x1
end;

plot(p->NRIG_quan(p,1),.01,.99)

@time NRIG_quan.(rand(10000),1);

  0.622641 seconds (2.48 M allocations: 60.534 MiB, 2.65% gc time)

Student's-t₄ random variable¶

"conventional Student's-t4 density"
function t4_dens(x)
    return 3/8*(1+x^2/4)^(-5/2)
end;

"conventional Student's-t₄ distribution"
function t4_dist(x)
    return (1+(x*(6+x^2))/((4+x^2)^(3/2)))/2
end;

plot([t4_dens,t4_dist],-2,2)

"conventional Student's-t₄ quantile"
function t4_quan(p)
    a=1/(4p*(1-p))
    return 2sign(p-1/2)sqrt(sqrt(a)cos(atan(sqrt(a-1))/3)-1)
end;

plot(t4_quan,.01,.99)

Normal random variable¶

symmetric Normal / reciprocal inverse Gaussian (NRIG) random variable¶

Student's-t4 random variable¶

Student's-t₄ random variable¶