diff --git a/gen_imgs.py b/gen_imgs.py
index 4e9d3ca0..6fc8694d 100644
--- a/gen_imgs.py
+++ b/gen_imgs.py
@@ -323,7 +323,6 @@ def chicago_tensor():
 
 
 def expectiles():
-
     """Generate expectiles visualization."""
     X, y = mcycle(return_X_y=True)
 
diff --git a/pygam/links.py b/pygam/links.py
index cc976eec..328fe726 100644
--- a/pygam/links.py
+++ b/pygam/links.py
@@ -118,7 +118,7 @@ def mu(self, lp, dist):
         -------
         mu : np.array of length n
         """
-        elp = np.exp(lp)
+        elp = np.exp(np.clip(lp, -500, 500))
         return dist.levels * elp / (elp + 1)
 
     def gradient(self, mu, dist):
@@ -178,7 +178,7 @@ def mu(self, lp, dist):
         -------
         mu : np.array of length n
         """
-        return np.exp(lp)
+        return np.exp(np.clip(lp, -500, 500))
 
     def gradient(self, mu, dist):
         """
diff --git a/pygam/tests/test_links.py b/pygam/tests/test_links.py
new file mode 100644
index 00000000..ae5fc40f
--- /dev/null
+++ b/pygam/tests/test_links.py
@@ -0,0 +1,21 @@
+import numpy as np
+
+from pygam.distributions import BinomialDist
+from pygam.links import LogitLink, LogLink
+
+
+def test_log_link_no_overflow():
+    """LogLink.mu should not overflow for large linear predictor values."""
+    link = LogLink()
+    lp = np.array([-1000.0, 0.0, 1000.0])
+    result = link.mu(lp, dist=None)
+    assert np.all(np.isfinite(result)), "LogLink.mu produced inf or nan"
+
+
+def test_logit_link_no_overflow():
+    """LogitLink.mu should not overflow for large linear predictor values."""
+    link = LogitLink()
+    dist = BinomialDist()
+    lp = np.array([-1000.0, 0.0, 1000.0])
+    result = link.mu(lp, dist)
+    assert np.all(np.isfinite(result)), "LogitLink.mu produced inf or nan"
diff --git a/test_output.txt b/test_output.txt
new file mode 100644
index 00000000..f2aa924a
Binary files /dev/null and b/test_output.txt differ