diff --git a/src/Utilities.jl b/src/Utilities.jl
index 563971b..7036102 100644
--- a/src/Utilities.jl
+++ b/src/Utilities.jl
@@ -84,18 +84,64 @@
         v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
         return exp(d.A + d.shp*d.skw*xs*v - d.shp/d.skw*xs*(1-v))
     end
-    Distributions.cdf(d::BilinearExponential, x::Real) = first(quadgk(x->pdf(d,x), d.μ-200*d.σ/d.shp, x))
-    Distributions.ccdf(d::BilinearExponential, x::Real) = first(quadgk(x->pdf(d,x), x, d.μ+200*d.σ/d.shp))
+    function Distributions.cdf(d::BilinearExponential{T}, x::Real) where {T}
+        maxspan = 200*d.σ/d.shp
+        if d.μ + maxspan < x
+            one(float(T))
+        elseif x < d.μ - maxspan
+            zero(float(T))
+        else
+            l = min(d.μ, x)
+            return first(quadgk(x->pdf(d,x), l-maxspan, x))
+        end
+    end
+    function Distributions.ccdf(d::BilinearExponential{T}, x::Real) where {T}
+        maxspan =  200*d.σ/d.shp
+        if x < d.μ - maxspan
+            one(float(T))
+        elseif d.μ + maxspan < x
+            zero(float(T))
+        else
+            u = max(d.μ, x)
+            first(quadgk(x->pdf(d,x), x, u+maxspan))
+        end
+    end
+
     @inline function Distributions.logpdf(d::BilinearExponential, x::Real)
         xs = (x - d.μ)/d.σ # X scaled by mean and variance
         v = 1/2 - atan(xs)/3.141592653589793 # Sigmoid (positive on LHS)
         return d.A + d.shp*d.skw*xs*v - d.shp/d.skw*xs*(1-v)
     end
+    function Distributions.logcdf(d::BilinearExponential{T}, x::Real) where {T}
+        maxspan = 200*d.σ/d.shp
+        if x > d.μ + maxspan
+            zero(float(T))
+        else
+            l = min(d.μ, x)
+            c = logpdf(d, l)
+            area_shifted, e = quadgk(x->exp(logpdf(d,x)-c), l-maxspan, x)
+            log(area_shifted)+c
+        end
+    end
+    function Distributions.logccdf(d::BilinearExponential{T}, x::Real) where {T}
+        maxspan = 200*d.σ/d.shp
+        if x < d.μ - maxspan
+            zero(float(T))
+        else
+            u = max(d.μ, x)
+            c = logpdf(d, u)
+            area_shifted, e = quadgk(x->exp(logpdf(d,x)-c), x, u+maxspan)
+            log(area_shifted)+c
+        end
+    end
+
 
     ## Statistics
     Distributions.mean(d::BilinearExponential) = first(quadgk(x->x*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))
     Distributions.var(d::BilinearExponential; mean=mean(d)) = first(quadgk(x->(x-mean)^2*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))
     Distributions.std(d::BilinearExponential; mean=mean(d)) = sqrt(var(d; mean))
+    Distributions.skewness(d::BilinearExponential; mean=mean(d), std=std(d; mean)) = first(quadgk(x->(x-mean)^3*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))/std^3
+    Distributions.kurtosis(d::BilinearExponential; mean=mean(d), std=std(d; mean)) = first(quadgk(x->(x-mean)^4*pdf(d,x), d.μ-100*d.σ/d.shp, d.μ+100*d.σ/d.shp, maxevals=1000))/std^4
 
     ## Affine transformations
     Base.:+(d::BilinearExponential{T}, c::Real) where {T} = BilinearExponential{T}(d.A, d.μ + c, d.σ, d.shp, d.skw)
diff --git a/test/testUtilities.jl b/test/testUtilities.jl
index f90d884..e0c28fb 100644
--- a/test/testUtilities.jl
+++ b/test/testUtilities.jl
@@ -5,8 +5,11 @@
     @test pdf(d, 66.02) ≈ 0.0005437680417927236
     @test logpdf(d, 66.02) ≈ (-3.251727600957773 -4.265260194845627)
     @test pdf.(d, [66.02, 66.08, 66.15]) ≈ [0.0005437680417927236, 5.505685686251538e-13, 5.610687893459488e-28]
-    @test cdf.(d, [66.02, 66.08, 66.15]) ≈ [0.9999979646869678, 0.9999999999999903, 1.0000000000000446]
-    @test ccdf.(d, [66.02, 66.08, 66.15]) ≈ [2.035313044102181e-6, 1.2768666553428429e-15, 1.0225965286520566e-30]
+    @test cdf.(d, [10, 66.02, 66.08, 66.15, 100]) ≈ [0, 0.9999979646869678, 0.9999999999999903, 1.0000000000000446, 1]
+    @test ccdf.(d, [10, 66.02, 66.08, 66.15, 100]) ≈ [1, 2.035313044102181e-6, 1.2768666553428429e-15, 1.0225965286520566e-30, 0]
+    @test logcdf.(d, [10, 66.02, 66.08, 66.15, 100]) ≈  [-13366.73582762445, -2.035315134207849e-6, 6.394884621840902e-14, 3.4638958368304884e-14, 0.0]
+    @test logccdf.(d, [10, 66.02, 66.08, 66.15, 100]) ≈ [0.0, -13.10486092087404, -34.29436724356088, -69.05520778079227, -22437.23853845797]
+
 
     @test eltype(d) === partype(d) === Float64
     @test params(d) == (-4.265260194845627, 66.00606672179812, 0.17739474265630253, 70.57882331291309, 0.6017142541555505)
@@ -20,6 +23,8 @@
     @test mean(d) ≈ 65.93219453443876 atol = 1e-6
     @test var(d) ≈ 0.02080666596420308^2 atol = 1e-6
     @test std(d) ≈ 0.02080666596420308 atol = 1e-6
+    @test skewness(d) ≈ -0.16449904171203025 atol = 1e-6
+    @test kurtosis(d) ≈ 3.0865922917156188 atol = 1e-6
 
 ## --- Test Radiocarbon distribution