diff --git a/src/probabilit/distributions.py b/src/probabilit/distributions.py
index 91d6d44..f15341d 100644
--- a/src/probabilit/distributions.py
+++ b/src/probabilit/distributions.py
@@ -40,23 +40,30 @@ def Triangular(low, mode, high, low_perc=0.1, high_perc=0.9):
     Distribution("triang", loc=-2.2360679774997894, scale=14.472135954999578, c=0.5000000000000001)
     >>> Triangular(low=1, mode=5, high=9, low_perc=0.25, high_perc=0.75)
     Distribution("triang", loc=-8.656854249492383, scale=27.313708498984766, c=0.5)
+    >>> Triangular(low=1, mode=5, high=9, low_perc=0, high_perc=1)
+    Distribution("triang", loc=1, scale=8, c=0.5)
     """
     # A few comments on fitting can be found here:
     # https://docs.analytica.com/index.php/Triangular10_50_90
 
     if not (low < mode < high):
         raise ValueError(f"Must have {low=} < {mode=} < {high=}")
-    if not ((0 < low_perc <= 1.0) and (0 <= high_perc < 1.0)):
+    if not ((0 <= low_perc <= 1.0) and (0 <= high_perc <= 1.0)):
         raise ValueError("Percentiles must be between 0 and 1.")
 
-    # Optimize parameters
-    loc, scale, c = _fit_triangular_distribution(
-        low=low,
-        mode=mode,
-        high=high,
-        low_perc=low_perc,
-        high_perc=high_perc,
-    )
+    # No need to optimize if low and high are boundaries of distribution support
+    if np.isclose(low_perc, 0.0) and np.isclose(high_perc, 1.0):
+        loc, scale, c = low, high - low, (mode - low) / (high - low)
+
+    else:
+        # Optimize parameters
+        loc, scale, c = _fit_triangular_distribution(
+            low=low,
+            mode=mode,
+            high=high,
+            low_perc=low_perc,
+            high_perc=high_perc,
+        )
     return Distribution("triang", loc=loc, scale=scale, c=c)
 
 
diff --git a/tests/test_modeling.py b/tests/test_modeling.py
index 5bd6d37..aabd160 100644
--- a/tests/test_modeling.py
+++ b/tests/test_modeling.py
@@ -113,10 +113,10 @@ def test_total_person_hours(self):
 
         for i in range(num_rivets):
             total_person_hours += Triangular(
-                low=3.75, mode=4.25, high=5.5, low_perc=0.00001, high_perc=0.99999
+                low=3.75, mode=4.25, high=5.5, low_perc=0, high_perc=1.0
             )
 
-        num_samples = 10000
+        num_samples = 1000
         res_total_person_hours = total_person_hours.sample(num_samples, rng)
 
         # The mean and standard deviation of a Triangular(3.75, 4.25, 5.5) are 4.5 and 0.368,
@@ -126,10 +126,11 @@ def test_total_person_hours(self):
         expected_std = 0.368 * np.sqrt(num_rivets)
 
         sample_mean = np.mean(res_total_person_hours)
-        sample_std = np.std(res_total_person_hours)
+        sample_std = np.std(res_total_person_hours, ddof=1)
 
-        assert abs(sample_mean - expected_mean) < 0.3
-        assert abs(sample_std - expected_std) < 0.1
+        # Within 2% of theoretical values
+        np.testing.assert_allclose(sample_mean, expected_mean, rtol=0.02)
+        np.testing.assert_allclose(sample_std, expected_std, rtol=0.02)
 
 
 def test_copying():