47 lines
1.0 KiB
C
47 lines
1.0 KiB
C
|
#include "common.h"
|
||
|
|
#include <gsl/gsl_randist.h>
|
||
|
|
|
||
|
|
/* Create a sample of `n` points */
|
||
|
|
sample_t* sample_t_alloc(size_t n, struct par p) {
|
||
|
|
sample_t *x = malloc(sizeof(sample_t));
|
||
|
|
x->p = p;
|
||
|
|
x->data = gsl_matrix_alloc(n, 2);
|
||
|
|
return x;
|
||
|
|
}
|
||
|
|
|
||
|
|
/* Delete a sample */
|
||
|
|
void sample_t_free(sample_t *x) {
|
||
|
|
gsl_matrix_free(x->data);
|
||
|
|
free(x);
|
||
|
|
}
|
||
|
|
|
||
|
|
|
||
|
|
/* `generate_normal(r, n, p)` will create
|
||
|
|
* a sample of `n` points, distributed
|
||
|
|
* according to a bivariate gaussian distribution
|
||
|
|
* of parameters `p`.
|
||
|
|
*/
|
||
|
|
sample_t* generate_normal(
|
||
|
|
gsl_rng *r, size_t n, struct par *p) {
|
||
|
|
sample_t *s = sample_t_alloc(n, *p);
|
||
|
|
|
||
|
|
for (size_t i = 0; i < n; i++) {
|
||
|
|
/* Generate a vector (x,y) with
|
||
|
|
* a standard (μ = 0) bivariate
|
||
|
|
* gaussian PDF.
|
||
|
|
*/
|
||
|
|
double *x = gsl_matrix_ptr(s->data, i, 0);
|
||
|
|
double *y = gsl_matrix_ptr(s->data, i, 1);
|
||
|
|
gsl_ran_bivariate_gaussian(
|
||
|
|
r,
|
||
|
|
p->sigma_x, p->sigma_y, p->rho,
|
||
|
|
x, y);
|
||
|
|
|
||
|
|
/* Shift the vector to (x₀,y₀) */
|
||
|
|
*x += p->x0;
|
||
|
|
*y += p->y0;
|
||
|
|
}
|
||
|
|
|
||
|
|
return s;
|
||
|
|
}
|