47 static unsigned int pf_pdf_seed;
55 pf_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);
83 void pf_pdf_gaussian_free(pf_pdf_gaussian_t *pdf)
114 pf_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];
141 double pf_ran_gaussian(
double sigma)
147 do { r = drand48(); }
while (r==0.0);
149 do { r = drand48(); }
while (r==0.0);
152 }
while(w > 1.0 || w==0.0);
154 return(sigma * x2 * sqrt(-2.0*log(w)/w));
168 pf_pdf_discrete_t *pf_pdf_discrete_alloc(
int count,
double *probs)
170 pf_pdf_discrete_t *pdf;
172 pdf = calloc(1,
sizeof(pf_pdf_discrete_t));
174 pdf->prob_count = count;
175 pdf->probs = malloc(count *
sizeof(
double));
176 memcpy(pdf->probs, probs, count *
sizeof(
double));
179 pdf->rng = gsl_rng_alloc(gsl_rng_taus);
180 gsl_rng_set(pdf->rng, ++pf_pdf_seed);
183 pdf->ran = gsl_ran_discrete_preproc(count, probs);
190 void pf_pdf_discrete_free(pf_pdf_discrete_t *pdf)
192 gsl_ran_discrete_free(pdf->ran);
193 gsl_rng_free(pdf->rng);
201 double pf_pdf_discrete_value(pf_pdf_discrete_t *pdf,
int i)
203 return pdf->probs[i];
208 int pf_pdf_discrete_sample(pf_pdf_discrete_t *pdf)
212 i = gsl_ran_discrete(pdf->rng, pdf->ran);
213 assert(i >= 0 && i < pdf->prob_count);