47static unsigned int pf_pdf_seed;
55pf_pdf_gaussian_t *pf_pdf_gaussian_alloc(pf_vector_t *x, pf_matrix_t *cx)
58 pf_pdf_gaussian_t *pdf;
60 pdf = calloc(1,
sizeof(pf_pdf_gaussian_t));
68 pf_matrix_unitary(&pdf->cr, &cd, &pdf->cx);
69 pdf->cd.v[0] = sqrt(cd.m[0][0]);
70 pdf->cd.v[1] = sqrt(cd.m[1][1]);
71 pdf->cd.v[2] = sqrt(cd.m[2][2]);
76 srand48(++pf_pdf_seed);
83void pf_pdf_gaussian_free(pf_pdf_gaussian_t *pdf)
114pf_vector_t pf_pdf_gaussian_sample(pf_pdf_gaussian_t *pdf)
121 for (i = 0; i < 3; i++)
124 r.v[i] = pf_ran_gaussian(pdf->cd.v[i]);
127 for (i = 0; i < 3; i++)
129 x.v[i] = pdf->x.v[i];
130 for (j = 0; j < 3; j++)
131 x.v[i] += pdf->cr.m[i][j] * r.v[j];
141double pf_ran_gaussian(
double sigma)
148 do { r = drand48(); }
while (r==0.0);
149 double x1 = 2.0 * r - 1.0;
150 do { r = drand48(); }
while (r==0.0);
153 }
while(w > 1.0 || w==0.0);
155 return(sigma * x2 * sqrt(-2.0*log(w)/w));
169pf_pdf_discrete_t *pf_pdf_discrete_alloc(
int count,
double *probs)
171 pf_pdf_discrete_t *pdf;
173 pdf = calloc(1,
sizeof(pf_pdf_discrete_t));
175 pdf->prob_count = count;
176 pdf->probs = malloc(count *
sizeof(
double));
177 memcpy(pdf->probs, probs, count *
sizeof(
double));
180 pdf->rng = gsl_rng_alloc(gsl_rng_taus);
181 gsl_rng_set(pdf->rng, ++pf_pdf_seed);
184 pdf->ran = gsl_ran_discrete_preproc(count, probs);
191void pf_pdf_discrete_free(pf_pdf_discrete_t *pdf)
193 gsl_ran_discrete_free(pdf->ran);
194 gsl_rng_free(pdf->rng);
202double pf_pdf_discrete_value(pf_pdf_discrete_t *pdf,
int i)
204 return pdf->probs[i];
209int pf_pdf_discrete_sample(pf_pdf_discrete_t *pdf)
213 i = gsl_ran_discrete(pdf->rng, pdf->ran);
214 assert(i >= 0 && i < pdf->prob_count);