aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--gaussian-nd-random.c51
-rw-r--r--gaussian-nd-random.h19
2 files changed, 70 insertions, 0 deletions
diff --git a/gaussian-nd-random.c b/gaussian-nd-random.c
new file mode 100644
index 0000000..16abc94
--- /dev/null
+++ b/gaussian-nd-random.c
@@ -0,0 +1,51 @@
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_sf_gamma.h>
+#include "nd-random.h"
+#include "utils.h"
+
+#define INTEGRATION_EPSABS 1e-3
+#define INTEGRATION_EPSREL 1e-3
+
+unsigned int shifted_gaussian_random_vector
+(const gsl_rng* r, const gsl_vector* mean, double theta_max, double standard_deviation, gsl_vector* x)
+{
+ for (int i=0; i<x->size; i++)
+ gsl_vector_set(x, i, gsl_ran_gaussian(r, standard_deviation));
+ gsl_vector_add(x, mean);
+ gsl_blas_dscal(1/gsl_blas_dnrm2(x), x);
+ if (angle_between_vectors(x, mean) > theta_max)
+ return 1 + shifted_gaussian_random_vector(r, mean, theta_max, standard_deviation, x);
+ else return 1;
+}
+
+double shifted_gaussian_pdf (double theta, double mean, double theta_max,
+ double standard_deviation, unsigned int dimension, gsl_integration_workspace* w)
+{
+ // NOTE: Due to the rejection sampling, the returned result is only
+ // proportional to the actual PDF.
+ if (theta > theta_max) return 0;
+
+ double projection = mean*cos(theta);
+ double integrand (double r, void* params)
+ {
+ return gsl_pow_uint(r, dimension - 1)
+ * gaussian_pdf((r - projection) / standard_deviation);
+ }
+
+ gsl_function gsl_integrand = {&integrand};
+ double result, error;
+ /* gsl_integration_qagiu(&gsl_integrand, 0, INTEGRATION_EPSABS, INTEGRATION_EPSREL, */
+ /* w->limit, w, &result, &error); */
+ double rmax = 0.5 * (projection + sqrt(gsl_pow_2(projection) + 4*(dimension-1)*gsl_pow_2(standard_deviation)));
+ double halfwidth = 3 * standard_deviation;
+ gsl_integration_qag(&gsl_integrand,
+ rmax - halfwidth < 0 ? 0 : rmax - halfwidth,
+ rmax + halfwidth,
+ INTEGRATION_EPSABS, INTEGRATION_EPSREL,
+ w->limit, GSL_INTEG_GAUSS41, w, &result, &error);
+ return result
+ * gaussian_pdf(sqrt(gsl_pow_2(mean) - gsl_pow_2(projection)) / standard_deviation)
+ / gsl_pow_uint(M_SQRT2*sqrt(M_PI), dimension - 2)
+ / gsl_pow_uint(standard_deviation, dimension);
+}
diff --git a/gaussian-nd-random.h b/gaussian-nd-random.h
new file mode 100644
index 0000000..951cb10
--- /dev/null
+++ b/gaussian-nd-random.h
@@ -0,0 +1,19 @@
+#ifndef GAUSSIAN_ND_RANDOM_H
+#define GAUSSIAN_ND_RANDOM_H
+
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_integration.h>
+
+double planar_angle_to_standard_deviation
+(double mean, double theta_max, double truncation, unsigned int dimension);
+
+unsigned int shifted_gaussian_random_vector
+(const gsl_rng* r, const gsl_vector* mean,
+ double theta_max, double truncation, gsl_vector* x);
+
+double shifted_gaussian_pdf
+(double theta, double mean, double theta_max,
+ double truncation, unsigned int dimension, gsl_integration_workspace* w);
+
+#endif