Support for zero σ normals (#62)

This adds support for normal distributions with zero standard deviation. At the moment, this gives the same values as R. The only one I'm not particularly sure about is what to return for the quantile at 1 (R gives +Inf, there is a reasonable argument I think it should be μ).

src/distrs/norm.jl

@@ -1,35 +1,92 @@
 # functions related to normal distribution
 
-xval(μ::Real, σ::Real, z::Number) = μ + σ * z
+function xval(μ::Real, σ::Real, z::Number)
+    if isinf(z) && iszero(σ)
+        μ + one(σ) * z
+    else
+        μ + σ * z
+    end
+end
 zval(μ::Real, σ::Real, x::Number) = (x - μ) / σ
 
 # pdf
 normpdf(z::Number) = exp(-abs2(z)/2) * invsqrt2π
-normpdf(μ::Real, σ::Real, x::Number) = normpdf(zval(μ, σ, x)) / σ
+function normpdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ)
+        if x == μ
+            z = zval(μ, one(σ), x)
+        else
+            z = zval(μ, σ, x)
+            σ = one(σ)
+        end
+    else
+        z = zval(μ, σ, x)
+    end
+    normpdf(z) / σ
+end
 
 # logpdf
 normlogpdf(z::Number) = -(abs2(z) + log2π)/2
-normlogpdf(μ::Real, σ::Real, x::Number) = normlogpdf(zval(μ, σ, x)) - log(σ)
+function normlogpdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ)
+        if x == μ
+            z = zval(μ, one(σ), x)
+        else
+            z = zval(μ, σ, x)
+            σ = one(σ)
+        end
+    else
+        z = zval(μ, σ, x)
+    end
+    normlogpdf(z) - log(σ)
+end            
 
 # cdf
 normcdf(z::Number) = erfc(-z * invsqrt2)/2
-normcdf(μ::Real, σ::Real, x::Number) = normcdf(zval(μ, σ, x))
-
+function normcdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ) && x == μ
+        z = zval(zero(μ), σ, one(x))
+    else        
+        z = zval(μ, σ, x)        
+    end
+    normcdf(z)
+end
 # ccdf
 normccdf(z::Number) = erfc(z * invsqrt2)/2
-normccdf(μ::Real, σ::Real, x::Number) = normccdf(zval(μ, σ, x))
+function normccdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ) && x == μ
+        z = zval(zero(μ), σ, one(x))
+    else        
+        z = zval(μ, σ, x)        
+    end
+    normccdf(z)
+end
 
 # logcdf
 normlogcdf(z::Number) = z < -1.0 ?
     log(erfcx(-z * invsqrt2)/2) - abs2(z)/2 :
     log1p(-erfc(z * invsqrt2)/2)
-normlogcdf(μ::Real, σ::Real, x::Number) = normlogcdf(zval(μ, σ, x))
+function normlogcdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ) && x == μ
+        z = zval(zero(μ), σ, one(x))
+    else        
+        z = zval(μ, σ, x)        
+    end
+    normlogcdf(z)
+end
 
 # logccdf
 normlogccdf(z::Number) = z > 1.0 ?
     log(erfcx(z * invsqrt2)/2) - abs2(z)/2 :
     log1p(-erfc(-z * invsqrt2)/2)
-normlogccdf(μ::Real, σ::Real, x::Number) = normlogccdf(zval(μ, σ, x))
+function normlogccdf(μ::Real, σ::Real, x::Number)
+    if iszero(σ) && x == μ
+        z = zval(zero(μ), σ, one(x))
+    else        
+        z = zval(μ, σ, x)        
+    end
+    normlogccdf(z)
+end
 
 norminvcdf(p::Real) = -erfcinv(2*p) * sqrt2
 norminvcdf(μ::Real, σ::Real, p::Real) = xval(μ, σ, norminvcdf(p))

test/rmath.jl

@@ -5,20 +5,9 @@ function check_rmath(fname, statsfun, rmathfun, params, aname, a, isprob, rtol)
     v = statsfun(params..., a)
     rv = rmathfun(params..., a)
     if isprob
-        rd = abs(v / rv - 1.0)
-        if rd > rtol
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v/rv - 1| = $rd > $rtol.")
-        end
+        @test v ≈ rv rtol=rtol nans=true
     else
-        τ = (1.0 + abs(rv)) * rtol
-        ad = abs(v - rv)
-        if ad > τ
-            error("$fname deviates too much from Rmath at " *
-                "params = $params, $aname = $a:\n" *
-                "  v = $v (rv = $rv)\n  |v - rv| = $ad > $τ.")
-        end
+        @test v ≈ rv atol=rtol rtol=rtol nans=true
     end
 end
 
@@ -67,8 +56,7 @@ function rmathcomp(basename, params, X::AbstractArray, rtol=100eps(float(one(elt
     end
     rmath_rand = has_rand ? get_rmathfun(rand) : nothing
 
-    for i = 1:length(X)
-        x = X[i]
+    for x in X
         check_rmath(pdf, stats_pdf, rmath_pdf,
             params, "x", x, true, rtol)
         check_rmath(logpdf, stats_logpdf, rmath_logpdf,
@@ -195,7 +183,8 @@ rmathcomp_tests("norm", [
     ((0.0, 0.5), -3.0:0.01:3.0),
     ((0, 1), -3.0:3.0),
     ((0, 1), -3.0:0.01:3.0),
-    ((0, 1), -3:3)
+    ((0, 1), -3:3),
+    ((0.0, 0.0), -6.0:0.1:6.0),
 ])
 
 rmathcomp_tests("ntdist", [