Use truncnorm for bounded Gaussian prior

CHANGELOG-unreleased.md

@@ -24,4 +24,5 @@ the released changes.
 - Make `get_prefix_timeranges` work for SWX.
 - Some of the `gridutils` functions had improper logging behavior
 - Fixed bug in changing epoch for ELL1k model
+- Fixed bug in `GaussianRV_gen`, where the probability distribution function was not normalized correctly. Changed to use `scipy.stats.truncnorm` instead of the custom `GaussianRV_gen`.
 ### Removed

src/pint/models/priors.py

@@ -8,7 +8,7 @@
 """
 
 import numpy as np
-from scipy.stats import rv_continuous, uniform
+from scipy.stats import rv_continuous, truncnorm, uniform
 
 
 class Prior:
@@ -145,10 +145,9 @@ def GaussianBoundedRV(loc=0.0, scale=1.0, lower_bound=-np.inf, upper_bound=np.in
     upper_bound : number
         Upper bound of allowed parameter range
 
-    Returns a frozen rv_continuous instance with a gaussian probability
-    inside the range lower_bound to upper_bound and 0.0 outside
+    Returns a frozen truncnorm instance with a gaussian probability inside
+    the range lower_bound to upper_bound and 0.0 outside
     """
     ymin = (lower_bound - loc) / scale
     ymax = (upper_bound - loc) / scale
-    n = GaussianRV_gen(name="bounded_gaussian", a=ymin, b=ymax)
-    return n(loc=loc, scale=scale)
+    return truncnorm(a=ymin, b=ymax, loc=loc, scale=scale)

tests/test_priors.py

@@ -72,15 +72,17 @@ def test_inclination_prior(self):
     def test_gaussian_bounded(self):
         print("test_gaussian_bounded")
         self.m.M2.prior = Prior(
-            GaussianBoundedRV(loc=0.26, scale=0.10, lower_bound=0.0, upper_bound=0.6)
+            GaussianBoundedRV(loc=0.26, scale=0.20, lower_bound=0.06, upper_bound=0.46)
         )
         assert self.m.M2.prior_pdf(-0.1) == 0.0
         assert self.m.M2.prior_pdf(0.7) == 0.0
         assert self.m.M2.prior_pdf(-0.1, logpdf=True) == -np.inf
         assert self.m.M2.prior_pdf(0.7, logpdf=True) == -np.inf
         assert np.isclose(
-            self.m.M2.prior_pdf(0.26 + 0.1) / self.m.M2.prior_pdf(0.26),
+            self.m.M2.prior_pdf(0.26 + 0.2) / self.m.M2.prior_pdf(0.26),
             0.60653065971263342,
         )
         # Test that integral is 1.0, not safe since using _rv private var
         assert np.isclose(self.m.M2.prior._rv.cdf(0.6), 1.0)
+        # Test that integral is 0.5 at the mean, not safe since using _rv private var
+        assert np.isclose(self.m.M2.prior._rv.cdf(0.26), 0.5)