from msunet =at= shellblade =dot= net

Please find attached an implementation of a truncated normal distribution
based on John Burkardt's presentation linked below.

Please note that I am neither a mathematician nor an expert in numerical
computing, so any feedback is appreciated. This is also my first contribution
to GSL (and a GNU project in general), so please excuse any inefficiencies.

Thank you,

Marc



Adds a truncated normal distribution as described in John Burkardt's
"The Truncated Normal Distribution":

https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf
---
 cdf/gauss.c           | 22 ++++++++++++++++
 cdf/gaussinv.c        | 20 ++++++++++++++
 cdf/gsl_cdf.h         |  6 +++++
 cdf/test.c            | 61 +++++++++++++++++++++++++++++++++++++++++++
 filter/Makefile.am    |  2 +-
 movstat/Makefile.am   |  2 +-
 randist/gauss.c       | 19 ++++++++++++++
 randist/gsl_randist.h |  3 +++
 randist/test.c        | 20 ++++++++++++++
 rstat/Makefile.am     |  2 +-
 10 files changed, 154 insertions(+), 3 deletions(-)

diff --git a/cdf/gauss.c b/cdf/gauss.c
index cdc8b0a3..110e58f1 100644
--- a/cdf/gauss.c
+++ b/cdf/gauss.c
@@ -349,3 +349,25 @@ gsl_cdf_gaussian_Q (const double x, const double sigma)
 {
   return gsl_cdf_ugaussian_Q (x / sigma);
 }
+
+double
+gsl_cdf_tgaussian_P (const double x, const double a, const double b, const double sigma)
+{
+  if (x <= a)
+    {
+      return 0.0;
+    }
+  if (x >= b)
+    {
+      return 1.0;
+    }
+  const double n = gsl_cdf_gaussian_P (x, sigma) - gsl_cdf_gaussian_P (a, sigma);
+  const double d = gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma);
+  return n / d;
+}
+
+double
+gsl_cdf_tgaussian_Q (const double x, const double a, const double b, const double sigma)
+{
+  return 1.0 - gsl_cdf_tgaussian_P (x, a, b, sigma);
+}
diff --git a/cdf/gaussinv.c b/cdf/gaussinv.c
index dcbdca1d..28c9866c 100644
--- a/cdf/gaussinv.c
+++ b/cdf/gaussinv.c
@@ -200,3 +200,23 @@ gsl_cdf_gaussian_Qinv (const double Q, const double sigma)
 {
   return sigma * gsl_cdf_ugaussian_Qinv (Q);
 }
+
+double
+gsl_cdf_tgaussian_Pinv (const double P, const double a, const double b, const double sigma)
+{
+  const double k
+    = gsl_cdf_gaussian_P (a, sigma)
+    + P * (gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma));
+
+  return gsl_cdf_gaussian_Pinv (k, sigma);
+}
+
+double
+gsl_cdf_tgaussian_Qinv (const double Q, const double a, const double b, const double sigma)
+{
+  const double k
+    = gsl_cdf_gaussian_P (b, sigma)
+    - Q * (gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma));
+
+  return gsl_cdf_gaussian_Pinv (k, sigma);
+}
diff --git a/cdf/gsl_cdf.h b/cdf/gsl_cdf.h
index 2bc3fed5..61379a82 100644
--- a/cdf/gsl_cdf.h
+++ b/cdf/gsl_cdf.h
@@ -46,6 +46,12 @@ double gsl_cdf_gaussian_Q (const double x, const double sigma);
 double gsl_cdf_gaussian_Pinv (const double P, const double sigma);
 double gsl_cdf_gaussian_Qinv (const double Q, const double sigma);
 
+double gsl_cdf_tgaussian_P (const double x, const double a, const double b, const double sigma);
+double gsl_cdf_tgaussian_Q (const double x, const double a, const double b, const double sigma);
+
+double gsl_cdf_tgaussian_Pinv (const double P, const double a, const double b, const double sigma);
+double gsl_cdf_tgaussian_Qinv (const double Q, const double a, const double b, const double sigma);
+
 double gsl_cdf_gamma_P (const double x, const double a, const double b);
 double gsl_cdf_gamma_Q (const double x, const double a, const double b);
 
diff --git a/cdf/test.c b/cdf/test.c
index 1b4cf777..65eeaa6e 100644
--- a/cdf/test.c
+++ b/cdf/test.c
@@ -39,6 +39,8 @@
 
 void test_ugaussian (void);
 void test_ugaussianinv (void);
+void test_tgaussian (void);
+void test_tgaussianinv (void);
 void test_exponential (void);
 void test_exponentialinv (void);
 void test_exppow (void);
@@ -398,6 +400,7 @@ main (void)
      Function values computed with PARI, 28 digits precision */
  
   test_ugaussian ();
+  test_tgaussian ();
   test_exponential ();
   test_exppow ();
   test_tdist ();
@@ -407,6 +410,7 @@ main (void)
   test_beta ();
 
   test_ugaussianinv ();
+  test_tgaussianinv ();
   test_exponentialinv ();
   test_gammainv ();
   test_chisqinv ();
@@ -535,6 +539,63 @@ void test_ugaussianinv (void) {
 }
 
 
+  /* Test values from John Burkardt, The Truncated Normal Distribution
+     Page 23. */
+
+void test_tgaussian (void)
+{
+  TEST (gsl_cdf_tgaussian_P, (-18.37, -50.0, 50.0, 25.0), 0.218418626765711416, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (37.962, -50.0, 50.0, 25.0), 0.956315769095676504, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (22.367, -50.0, 50.0, 25.0), 0.82951388212098065, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (3.704, -50.0, 50.0, 25.0), 0.561699074822187949, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (-5.1, -50.0, 50.0, 25.0), 0.415323967636885893, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (-34.1674, -50.0, 50.0, 25.0), 0.0661185873044326383, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (-15.4257, -50.0, 50.0, 25.0), 0.257577857199668192, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (-28.4328, -50.0, 50.0, 25.0), 0.109956874884255498, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (-37.9346, -50.0, 50.0, 25.0), 0.0438289787355652244, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_P, (8.155, -50.0, 50.0, 25.0), 0.633958632430714042, TEST_TOL0);
+
+  TEST (gsl_cdf_tgaussian_Q, (-18.37, -50.0, 50.0, 25.0), 0.7815813732342887, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (37.962, -50.0, 50.0, 25.0), 0.043684230904323385, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (22.367, -50.0, 50.0, 25.0), 0.17048611787901935, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (3.704, -50.0, 50.0, 25.0), 0.43830092517781205, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (-5.1, -50.0, 50.0, 25.0), 0.584676032363114162, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (-34.1674, -50.0, 50.0, 25.0), 0.9338814126955673, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (-15.4257, -50.0, 50.0, 25.0), 0.7424221428003318, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (-28.4328, -50.0, 50.0, 25.0), 0.8900431251157445, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (-37.9346, -50.0, 50.0, 25.0), 0.9561710212644348, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Q, (8.155, -50.0, 50.0, 25.0), 0.36604136756928596, TEST_TOL0);
+}
+
+  /* Test values from John Burkardt, The Truncated Normal Distribution
+     Page 24. */
+
+void test_tgaussianinv (void)
+{
+  TEST (gsl_cdf_tgaussian_Pinv, (0.218418626765711305, -50.0, 50.0, 25.0), -18.37, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.956315769095676504, -50.0, 50.0, 25.0), 37.961999999999982, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.82951388212098065, -50.0, 50.0, 25.0), 22.367, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.561699074822187949, -50.0, 50.0, 25.0), 3.7039999999999913, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.415307593603450431, -50.0, 50.0, 25.0), -5.1009999999999982, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.0661185873044326383, -50.0, 50.0, 25.0), -34.1674, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.257577857199668192, -50.0, 50.0, 25.0), -15.4257, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.109956874884255498, -50.0, 50.0, 25.0), -28.4328, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.0438289787355652244, -50.0, 50.0, 25.0), -37.9346, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Pinv, (0.633958632430714042, -50.0, 50.0, 25.0), 8.155, TEST_TOL0);
+
+  TEST (gsl_cdf_tgaussian_Qinv, (0.7815813732342887, -50.0, 50.0, 25.0), -18.37, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.043684230904323496, -50.0, 50.0, 25.0), 37.961999999999982, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.17048611787901935, -50.0, 50.0, 25.0), 22.3669999999999902, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.43830092517781205, -50.0, 50.0, 25.0), 3.7039999999999913, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.5846924063965495, -50.0, 50.0, 25.0), -5.10099999999999554, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.9338814126955673, -50.0, 50.0, 25.0), -34.1674, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.7424221428003318, -50.0, 50.0, 25.0), -15.4257000000000151, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.8900431251157445, -50.0, 50.0, 25.0), -28.4328, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.9561710212644348, -50.0, 50.0, 25.0), -37.9345999999999819, TEST_TOL0);
+  TEST (gsl_cdf_tgaussian_Qinv, (0.36604136756928596, -50.0, 50.0, 25.0), 8.155, TEST_TOL0);
+}
+
+
   /* Tests for exponential cumulative distribution function
      Function values computed with PARI, 28 digits precision */
 
diff --git a/filter/Makefile.am b/filter/Makefile.am
index 03684916..fa426c5d 100644
--- a/filter/Makefile.am
+++ b/filter/Makefile.am
@@ -12,4 +12,4 @@ check_PROGRAMS = test
 TESTS = $(check_PROGRAMS)
 
 test_SOURCES = test.c
-test_LDADD = libgslfilter.la ../movstat/libgslmovstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../vector/libgslvector.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../block/libgslblock.la ../poly/libgslpoly.la ../sys/libgslsys.la ../utils/libutils.la
+test_LDADD = libgslfilter.la ../movstat/libgslmovstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../vector/libgslvector.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../block/libgslblock.la ../poly/libgslpoly.la ../cdf/libgslcdf.la ../sys/libgslsys.la ../utils/libutils.la
diff --git a/movstat/Makefile.am b/movstat/Makefile.am
index 851b241b..3abd7ead 100644
--- a/movstat/Makefile.am
+++ b/movstat/Makefile.am
@@ -33,4 +33,4 @@ check_PROGRAMS = test
 TESTS = $(check_PROGRAMS)
 
 test_SOURCES = test.c
-test_LDADD = libgslmovstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../vector/libgslvector.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../block/libgslblock.la ../sys/libgslsys.la ../utils/libutils.la
+test_LDADD = libgslmovstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../vector/libgslvector.la ../blas/libgslblas.la ../cblas/libgslcblas.la ../block/libgslblock.la ../cdf/libgslcdf.la ../sys/libgslsys.la ../utils/libutils.la
diff --git a/randist/gauss.c b/randist/gauss.c
index 1b3442ea..1a072186 100644
--- a/randist/gauss.c
+++ b/randist/gauss.c
@@ -20,6 +20,7 @@
 
 #include <config.h>
 #include <math.h>
+#include <gsl/gsl_cdf.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_rng.h>
 #include <gsl/gsl_randist.h>
@@ -141,3 +142,21 @@ gsl_ran_ugaussian_pdf (const double x)
   return gsl_ran_gaussian_pdf (x, 1.0);
 }
 
+double
+gsl_ran_tgaussian (const gsl_rng * r, const double a, const double b, const double sigma)
+{
+  const double p = gsl_rng_uniform_pos (r);
+  return gsl_cdf_tgaussian_Pinv (p, a, b, sigma);
+}
+
+double
+gsl_ran_tgaussian_pdf (const double x, const double a, const double b, const double sigma)
+{
+  if (x <= a || x >= b)
+    {
+      return 0.0;
+    }
+  const double n = gsl_ran_gaussian_pdf (x, sigma);
+  const double d = gsl_cdf_gaussian_P (b, sigma) - gsl_cdf_gaussian_P (a, sigma);
+  return n / d;
+}
diff --git a/randist/gsl_randist.h b/randist/gsl_randist.h
index d38ccb36..69027212 100644
--- a/randist/gsl_randist.h
+++ b/randist/gsl_randist.h
@@ -92,6 +92,9 @@ double gsl_ran_gaussian_tail_pdf (const double x, const double a, const double s
 double gsl_ran_ugaussian_tail (const gsl_rng * r, const double a);
 double gsl_ran_ugaussian_tail_pdf (const double x, const double a);
 
+double gsl_ran_tgaussian (const gsl_rng * r, const double a, const double b, const double sigma);
+double gsl_ran_tgaussian_pdf (const double x, const double a, const double b, const double sigma);
+
 void gsl_ran_bivariate_gaussian (const gsl_rng * r, double sigma_x, double sigma_y, double rho, double *x, double *y);
 double gsl_ran_bivariate_gaussian_pdf (const double x, const double y, const double sigma_x, const double sigma_y, const double rho);
 
diff --git a/randist/test.c b/randist/test.c
index c841f8b9..1b644259 100644
--- a/randist/test.c
+++ b/randist/test.c
@@ -156,6 +156,8 @@ double test_ugaussian_ratio_method (void);
 double test_ugaussian_ratio_method_pdf (double x);
 double test_ugaussian_tail (void);
 double test_ugaussian_tail_pdf (double x);
+double test_tgaussian (void);
+double test_tgaussian_pdf (double x);
 double test_bivariate_gaussian1 (void);
 double test_bivariate_gaussian1_pdf (double x);
 double test_bivariate_gaussian2 (void);
@@ -282,6 +284,11 @@ main (void)
   testMoments (FUNC (ugaussian), -1.0, 1.0, 0.6826895);
   testMoments (FUNC (ugaussian), 3.0, 3.5, 0.0011172689);
   testMoments (FUNC (ugaussian_tail), 3.0, 3.5, 0.0011172689 / 0.0013498981);
+  testMoments (FUNC (tgaussian), -3.0, 3.0, 1.0);
+  testMoments (FUNC (tgaussian), -9.0, -3.1, 0.0);
+  testMoments (FUNC (tgaussian), 3.1, 9.0, 0.0);
+  testMoments (FUNC (tgaussian), 0.0, 3.0, 0.5);
+  testMoments (FUNC (tgaussian), -3.0, 0.0, 0.5);
   testMoments (FUNC (exponential), 0.0, 1.0, 1 - exp (-0.5));
   testMoments (FUNC (cauchy), 0.0, 10000.0, 0.5);
 
@@ -339,6 +346,7 @@ main (void)
   testPDF (FUNC2 (gaussian_tail1));
   testPDF (FUNC2 (gaussian_tail2));
   testPDF (FUNC2 (ugaussian_tail));
+  testPDF (FUNC2 (tgaussian));
 
   testPDF (FUNC2 (bivariate_gaussian1));
   testPDF (FUNC2 (bivariate_gaussian2));
@@ -1626,6 +1634,18 @@ test_ugaussian_tail_pdf (double x)
   return gsl_ran_ugaussian_tail_pdf (x, 3.0);
 }
 
+double
+test_tgaussian (void)
+{
+  return gsl_ran_tgaussian (r_global, -3.0, 3.0, 1.5);
+}
+
+double
+test_tgaussian_pdf (double x)
+{
+  return gsl_ran_tgaussian_pdf (x, -3.0, 3.0, 1.5);
+}
+
 double
 test_bivariate_gaussian1 (void)
 {
diff --git a/rstat/Makefile.am b/rstat/Makefile.am
index 85db2013..5dd18e7b 100644
--- a/rstat/Makefile.am
+++ b/rstat/Makefile.am
@@ -10,6 +10,6 @@ check_PROGRAMS = test
 TESTS = $(check_PROGRAMS)
 
 test_SOURCES = test.c
-test_LDADD = libgslrstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la ../vector/libgslvector.la
+test_LDADD = libgslrstat.la ../statistics/libgslstatistics.la ../sort/libgslsort.la ../ieee-utils/libgslieeeutils.la ../randist/libgslrandist.la ../rng/libgslrng.la ../specfunc/libgslspecfunc.la ../complex/libgslcomplex.la ../err/libgslerr.la ../test/libgsltest.la ../sys/libgslsys.la ../utils/libutils.la ../vector/libgslvector.la ../cdf/libgslcdf.la
 
 
--
2.27.0