diff --git a/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py b/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py
index f4dcebb91..16e22ae3c 100644
--- a/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py
+++ b/src/pyrecest/distributions/hypersphere_subset/watson_distribution.py
@@ -155,7 +155,7 @@ def pdf(self, xs):
             raise ValueError(
                 f"xs must have trailing dimension {self.input_dim}, got {xs.shape}."
             )
-        p = self.norm_const * exp(self.kappa * (xs @ self.mu) ** 2)
+        p = exp(self.ln_norm_const + self.kappa * (xs @ self.mu) ** 2)
         return p
 
     def ln_pdf(self, xs):
diff --git a/tests/distributions/test_watson_distribution.py b/tests/distributions/test_watson_distribution.py
index 0b746ea84..004170bc7 100644
--- a/tests/distributions/test_watson_distribution.py
+++ b/tests/distributions/test_watson_distribution.py
@@ -81,6 +81,12 @@ def test_pdf(self):
         pdf_values = w.pdf(self.xs)
         npt.assert_array_almost_equal(pdf_values, expected_pdf_values, decimal=5)
 
+    def test_pdf_stays_finite_for_large_positive_kappa(self):
+        mu = array([1.0, 0.0, 0.0])
+        dist = WatsonDistribution(mu, 1000.0)
+
+        npt.assert_allclose(dist.pdf(mu), 159.07528584318962, rtol=1e-6)
+
     def test_pdf_accepts_list_inputs(self):
         mu = array([1.0, 2.0, 3.0])
         mu = mu / linalg.norm(mu)