star-sp

commit 197b852687fec5816e0eaa6370bcbdb93e15fef2
parent 66ce67f7e30b0f25b3ebc598af555bb8ab9b152b
Author: Christophe Coustet <christophe.coustet@meso-star.com>
Date:   Wed, 19 May 2021 12:56:53 +0200

Add truncated exponential distribution

Diffstat:
M
src/ssp.h
 | 
25
+++++++++++++++++++++++++
M
src/ssp_ran.c
 | 
44
++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 69 insertions(+), 0 deletions(-)
diff --git a/src/ssp.h b/src/ssp.h
@@ -327,6 +327,31 @@ ssp_ran_exp_float_pdf
   (const float x,
    const float mu);
 
+/* Truncated exponential distribution */
+SSP_API double
+ssp_ran_exp_truncated
+  (struct ssp_rng* rng,
+   const double mu,
+   const float max);
+
+SSP_API float
+ssp_ran_exp_truncated_float
+  (struct ssp_rng* rng,
+   const float mu,
+   const float max);
+
+SSP_API double
+ssp_ran_exp_truncated_pdf
+  (const double x,
+   const double mu,
+   const double max);
+
+SSP_API float
+ssp_ran_exp_truncated_float_pdf
+  (const float x,
+   const float mu,
+   const float max);
+
 /* Stateless random variate from the gaussian (or normal) distribution
  * with mean `mu':
  * P(x) dx = 1 / (sigma*sqrt(2*PI)) * exp(-1/2*((x-mu)/sigma)^2) dx */
diff --git a/src/ssp_ran.c b/src/ssp_ran.c
@@ -79,6 +79,50 @@ ssp_ran_exp_float_pdf(const float x, const float mu)
 }
 
 double
+ssp_ran_exp_truncated(struct ssp_rng* rng, const double mu, const double max)
+{
+  double u, norm;
+  ASSERT(rng && mu > 0 && max > 0);
+  u = ssp_rng_canonical(rng);
+  norm = 1 - exp(-mu * max);
+  return -log(norm * u) / mu;
+}
+
+float
+ssp_ran_exp_truncated_float(struct ssp_rng* rng, const float mu, const float max)
+{
+  float u, norm;
+  ASSERT(rng && mu > 0 && max > 0);
+  u = ssp_rng_canonical_float(rng);
+  norm = 1 - expf(-mu * max);
+  return -logf(norm * u) / mu;
+}
+
+double
+ssp_ran_exp_truncated_pdf(const double x, const double mu, const double max)
+{
+  ASSERT(x >= 0 && mu > 0 && max > 0);
+  if(x > max)
+    return 0;
+  else {
+    double norm = mu / (1 - exp(-max * mu));
+    return norm * exp(-x * mu);
+  }
+}
+
+float
+ssp_ran_exp_truncated_float_pdf(const float x, const float mu, const float max)
+{
+  ASSERT(x >= 0 && mu > 0 && max > 0);
+  if(x > max)
+    return 0;
+  else {
+    float norm = mu / (1 - expf(-max * mu));
+    return norm * expf(-x * mu);
+  }
+}
+
+double
 ssp_ran_gaussian
   (struct ssp_rng* rng,
    const double mu,