diff --git a/src/kaplanmeier.jl b/src/kaplanmeier.jl
index 5be057f..a63308a 100644
--- a/src/kaplanmeier.jl
+++ b/src/kaplanmeier.jl
@@ -41,9 +41,16 @@ function as a vector of tuples.
 function StatsAPI.confint(km::KaplanMeier; level::Real=0.05)
     q = quantile(Normal(), 1 - level/2)
     return map(km.survival, km.stderr) do srv, se
-        l = log(-log(srv))
+        # The direct implementation here would be
+        #   ℓ = log(-log(srv))
+        #   a = q * se / log(srv)
+        #   lower = exp(-exp(ℓ - a))
+        #   upper = exp(-exp(ℓ + a))
+        # However, this has some issues with numerical accuracy. An approximation to
+        # this quantity with improved accuracy is implemented here. The approximation
+        # was obtained via Herbie (https://herbie.uwplse.org/).
         a = q * se / log(srv)
-        exp(-exp(l - a)), exp(-exp(l + a))
+        return (srv^exp(-a), srv^exp(a))
     end
 end
 
diff --git a/src/nelsonaalen.jl b/src/nelsonaalen.jl
index 2c17676..2274ef7 100644
--- a/src/nelsonaalen.jl
+++ b/src/nelsonaalen.jl
@@ -40,8 +40,9 @@ function as a vector of tuples.
 """
 function StatsAPI.confint(na::NelsonAalen; level::Real=0.05)
     q = quantile(Normal(), 1 - level/2)
-    return map(na.chaz, na.stderr) do srv, se
-        srv - q * se, srv + q * se
+    return map(na.chaz, na.stderr) do ch, se
+        a = q * se
+        return (ch - a, ch + a)
     end
 end