1 /* Statistical routines for Weibull distributions.
2 *
3 * Contents:
4 * 1. Routines for evaluating densities and distributions.
5 * 2. Generic API routines, for general interface w/ histogram module
6 * 3. Dumping plots to files
7 * 4. Sampling
8 * 5. ML fitting to complete data
9 * 6. ML fitting to binned data
10 * 7. Test driver
11 * 8. Example
12 *
13 * To-do:
14 * - Fit*() functions should return eslEINVAL on n=0, eslENORESULT
15 * on failure due to small n. Compare esl_gumbel. xref J12/93.
16 * SRE, Wed Nov 27 11:13:48 2013
17 */
18 #include "esl_config.h"
19
20 #include <stdio.h>
21 #include <math.h>
22
23 #include "easel.h"
24 #include "esl_histogram.h"
25 #include "esl_minimizer.h"
26 #include "esl_random.h"
27 #include "esl_stats.h"
28 #include "esl_vectorops.h"
29
30 #include "esl_weibull.h"
31
32
33 /****************************************************************************
34 * 1. Routines for evaluating densities and distributions
35 ****************************************************************************/
36 /* mu <= x < infinity
37 * However, x=mu can be a problem:
38 * PDF-> 0 if tau > 1, infinity if tau < 1.
39 *
40 * lambda > 0
41 * tau > 0 [fat tail when tau < 1; inverse GEV when tau > 1;
42 * exponential when tau=1]
43 */
44
45
46 /* Function: esl_wei_pdf()
47 *
48 * Purpose: Calculates the Weibull pdf $P(X=x)$, given quantile <x>,
49 * offset <mu>, and parameters <lambda> and <tau>.
50 */
51 double
esl_wei_pdf(double x,double mu,double lambda,double tau)52 esl_wei_pdf(double x, double mu, double lambda, double tau)
53 {
54 double y = lambda * (x-mu);
55 double val;
56
57 if (x < mu) return 0.;
58 if (x == mu) {
59 if (tau < 1.) return eslINFINITY;
60 else if (tau > 1.) return 0.;
61 else if (tau == 1.) return lambda;
62 }
63
64 val = lambda * tau *
65 exp((tau-1)*log(y)) *
66 exp(- exp(tau * log(y)));
67 return val;
68 }
69
70 /* Function: esl_wei_logpdf()
71 *
72 * Purpose: Calculates the log probability density function for the
73 * Weibull, $\log P(X=x)$, given quantile <x>,
74 * offset <mu>, and parameters <lambda> and <tau>.
75 */
76 double
esl_wei_logpdf(double x,double mu,double lambda,double tau)77 esl_wei_logpdf(double x, double mu, double lambda, double tau)
78 {
79 double y = lambda * (x-mu);
80 double val;
81
82 if (x < mu) return -eslINFINITY;
83 if (x == mu) {
84 if (tau < 1.) return eslINFINITY; /* technically; but approaches it slowly*/
85 else if (tau > 1.) return -eslINFINITY; /* same as above, also a slow approach */
86 else if (tau == 1.) return log(lambda); /* special case, exponential */
87 }
88
89 val = log(tau) + tau*log(lambda) + (tau-1)*log(x-mu) - exp(tau * log(y));
90 return val;
91 }
92
93 /* Function: esl_wei_cdf()
94 *
95 * Purpose: Calculates the cumulative distribution function for the
96 * Weibull, $P(X \leq x)$, given quantile <x>,
97 * offset <mu>, and parameters <lambda> and <tau>.
98 */
99 double
esl_wei_cdf(double x,double mu,double lambda,double tau)100 esl_wei_cdf(double x, double mu, double lambda, double tau)
101 {
102 double y = lambda*(x-mu);
103 double tly = tau * log(y);
104
105 if (x <= mu) return 0.0;
106 else if (fabs(tly) < eslSMALLX1) return exp(tly);
107 else return 1 - exp(-exp(tly));
108 }
109
110 /* Function: esl_wei_logcdf()
111 *
112 * Purpose: Calculates the log of the cumulative distribution function for a
113 * Weibull, $P(X \leq x)$, given quantile <x>,
114 * offset <mu>, and parameters <lambda> and <tau>.
115 */
116 double
esl_wei_logcdf(double x,double mu,double lambda,double tau)117 esl_wei_logcdf(double x, double mu, double lambda, double tau)
118 {
119 double y = lambda*(x-mu);
120 double tly = tau * log(y);
121
122 if (x <= mu) return -eslINFINITY;
123
124 if (fabs(tly) < eslSMALLX1) return tly;
125 else if (fabs(exp(-exp(tly))) < eslSMALLX1) return -exp(-exp(tly));
126 else return log(1 - exp(-exp(tly)));
127 }
128
129
130 /* Function: esl_wei_surv()
131 *
132 * Purpose: Calculates the survivor function, $P(X>x)$ (that is, 1-CDF,
133 * the right tail probability mass) for a Weibull
134 * distribution, given quantile <x>, offset <mu>, and parameters
135 * <lambda> and <tau>.
136 */
137 double
esl_wei_surv(double x,double mu,double lambda,double tau)138 esl_wei_surv(double x, double mu, double lambda, double tau)
139 {
140 double y = lambda*(x-mu);
141 double tly = tau * log(y);
142
143 if (x <= mu) return 1.0;
144
145 return exp(-exp(tly));
146 }
147
148 /* Function: esl_wei_logsurv()
149 *
150 * Purpose: Calculates the log survivor function, $\log P(X>x)$ (that is,
151 * log(1-CDF), the right tail log probability mass) for a
152 * Weibull distribution, given quantile <x>, offset <mu>,
153 * and parameters <lambda> and <tau>.
154 */
155 double
esl_wei_logsurv(double x,double mu,double lambda,double tau)156 esl_wei_logsurv(double x, double mu, double lambda, double tau)
157 {
158 double y = lambda*(x-mu);
159 double tly = tau * log(y);
160
161 if (x <= mu) return 0.0;
162
163 return -exp(tly);
164 }
165
166 /* Function: esl_wei_invcdf()
167 *
168 * Purpose: Calculates the inverse CDF for a Weibull distribution
169 * with parameters <mu>, <lambda>, and <tau>, returning
170 * the quantile <x> at which the CDF is <p>, for $0<p<1$.
171 */
172 double
esl_wei_invcdf(double p,double mu,double lambda,double tau)173 esl_wei_invcdf(double p, double mu, double lambda, double tau)
174 {
175 return mu + 1/lambda * exp(1/tau * log(-log((1.-p))));
176 }
177 /*-------------------- end densities & distributions ------------------------*/
178
179
180
181
182 /****************************************************************************
183 * 2. Generic API routines: for general interface w/ histogram module
184 ****************************************************************************/
185
186 /* Function: esl_wei_generic_pdf()
187 *
188 * Purpose: Generic-API wrapper around <esl_wei_pdf()>, taking
189 * a void ptr to a double array containing $\mu$, $\lambda$,
190 * $\tau$ parameters.
191 */
192 double
esl_wei_generic_pdf(double x,void * params)193 esl_wei_generic_pdf(double x, void *params)
194 {
195 double *p = (double *) params;
196 return esl_wei_pdf(x, p[0], p[1], p[2]);
197 }
198
199 /* Function: esl_wei_generic_cdf()
200 *
201 * Purpose: Generic-API wrapper around <esl_wei_cdf()>, taking
202 * a void ptr to a double array containing $\mu$, $\lambda$,
203 * $\tau$ parameters.
204 */
205 double
esl_wei_generic_cdf(double x,void * params)206 esl_wei_generic_cdf(double x, void *params)
207 {
208 double *p = (double *) params;
209 return esl_wei_cdf(x, p[0], p[1], p[2]);
210 }
211
212 /* Function: esl_wei_generic_surv()
213 *
214 * Purpose: Generic-API wrapper around <esl_wei_surv()>, taking
215 * a void ptr to a double array containing $\mu$, $\lambda$,
216 * $\tau$ parameters.
217 */
218 double
esl_wei_generic_surv(double x,void * params)219 esl_wei_generic_surv(double x, void *params)
220 {
221 double *p = (double *) params;
222 return esl_wei_surv(x, p[0], p[1], p[2]);
223 }
224
225 /* Function: esl_wei_generic_invcdf()
226 *
227 * Purpose: Generic-API wrapper around <esl_wei_invcdf()>, taking
228 * a void ptr to a double array containing $\mu$, $\lambda$,
229 * $\tau$ parameters.
230 */
231 double
esl_wei_generic_invcdf(double p,void * params)232 esl_wei_generic_invcdf(double p, void *params)
233 {
234 double *v = (double *) params;
235 return esl_wei_invcdf(p, v[0], v[1], v[2]);
236 }
237 /*------------------------ end generic API ---------------------------------*/
238
239
240
241 /****************************************************************************
242 * 3. Dumping plots for files
243 ****************************************************************************/
244
245 /* Function: esl_wei_Plot()
246 *
247 * Purpose: Plot some Weibull function <func> (for instance, <esl_wei_pdf()>)
248 * for Weibull parameters <mu>, <lambda>, and <tau>, for a range of
249 * quantiles x from <xmin> to <xmax> in steps of <xstep>;
250 * output to an open stream <fp> in xmgrace XY input format.
251 *
252 * Returns: <eslOK> on success.
253 *
254 * Throws: <eslEWRITE> on any system write error, such as filled disk.
255 */
256 int
esl_wei_Plot(FILE * fp,double mu,double lambda,double tau,double (* func)(double x,double mu,double lambda,double tau),double xmin,double xmax,double xstep)257 esl_wei_Plot(FILE *fp, double mu, double lambda, double tau,
258 double (*func)(double x, double mu, double lambda, double tau),
259 double xmin, double xmax, double xstep)
260 {
261 double x;
262 for (x = xmin; x <= xmax; x += xstep)
263 if (x > mu || tau >= 1.) /* don't try to plot at mu where pdf blows up */
264 if (fprintf(fp, "%f\t%g\n", x, (*func)(x, mu, lambda, tau)) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "weibull plot write failed");
265 if (fprintf(fp, "&\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "weibull plot write failed");
266 return eslOK;
267 }
268 /*-------------------- end plot dumping routines ---------------------------*/
269
270
271
272
273
274 /****************************************************************************
275 * 4. Sampling
276 ****************************************************************************/
277
278 /* Function: esl_wei_Sample()
279 *
280 * Purpose: Sample a Weibull random variate,
281 * by the transformation method.
282 */
283 double
esl_wei_Sample(ESL_RANDOMNESS * r,double mu,double lambda,double tau)284 esl_wei_Sample(ESL_RANDOMNESS *r, double mu, double lambda, double tau)
285 {
286 double p;
287 p = esl_rnd_UniformPositive(r);
288 return esl_wei_invcdf(p, mu, lambda, tau);
289 }
290
291 /*--------------------------- end sampling ---------------------------------*/
292
293
294 /****************************************************************************
295 * 5. ML fitting to complete data
296 ****************************************************************************/
297
298 /* Easel's conjugate gradient descent code allows a single void ptr to
299 * point to any necessary fixed data, so we put everything into one
300 * structure:
301 */
302 struct wei_data {
303 double *x; /* data: n observed samples */
304 int n; /* number of observed samples */
305 double mu; /* mu is considered to be known, not fitted */
306 };
307
308 /* wei_func():
309 * Returns the negative log likelihood of a complete data sample,
310 * in the API of the conjugate gradient descent optimizer in esl_minimizer.
311 */
312 static double
wei_func(double * p,int nparam,void * dptr)313 wei_func(double *p, int nparam, void *dptr)
314 {
315 double lambda, tau;
316 struct wei_data *data;
317 double logL;
318 int i;
319
320 /* Unpack what the optimizer gave us.
321 */
322 lambda = exp(p[0]); /* see below for c.o.v. notes */
323 tau = exp(p[1]);
324 data = (struct wei_data *) dptr;
325
326 logL = 0.;
327 for (i = 0; i < data->n; i++)
328 {
329 if (tau < 1. && data->x[i] == data->mu) continue; /* hack: disallow infinity */
330 logL += esl_wei_logpdf(data->x[i], data->mu, lambda, tau);
331 }
332 return -logL; /* goal: minimize NLL */
333 }
334
335 /* Function: esl_wei_FitComplete()
336 *
337 * Purpose: Given an array of <n> samples <x[0]..x[n-1>, fit
338 * them to a stretched exponential distribution starting
339 * at lower bound <mu> (all $x_i > \mu$), and
340 * return maximum likelihood parameters <ret_lambda>
341 * and <ret_tau>.
342 *
343 * Args: x - complete GEV-distributed data [0..n-1]
344 * n - number of samples in <x>
345 * ret_mu - RETURN: lower bound of the distribution (all x_i >= mu)
346 * ret_lambda - RETURN: maximum likelihood estimate of lambda
347 * ret_tau - RETURN: maximum likelihood estimate of tau
348 *
349 * Returns: <eslOK> on success.
350 *
351 * Throws: <eslENOHALT> if the fit doesn't converge.
352 *
353 * Xref: STL9/136-137
354 */
355 int
esl_wei_FitComplete(double * x,int n,double * ret_mu,double * ret_lambda,double * ret_tau)356 esl_wei_FitComplete(double *x, int n, double *ret_mu,
357 double *ret_lambda, double *ret_tau)
358 {
359 struct wei_data data;
360 double p[2]; /* parameter vector */
361 double mean;
362 double mu, lambda, tau; /* initial param guesses */
363 double fx; /* f(x) at minimum; currently unused */
364 int status;
365
366 /* Make a good initial guess, based on exponential fit;
367 * set an arbitrary tau.
368 */
369 mu = esl_vec_DMin(x, n);
370 esl_stats_DMean(x, n, &mean, NULL);
371 lambda = 1 / (mean - mu);
372 tau = 0.9;
373
374 /* Load the data structure
375 */
376 data.x = x;
377 data.n = n;
378 data.mu = mu;
379
380 /* Change of variables;
381 * lambda > 0, so c.o.v. lambda = exp^w, w = log(lambda);
382 * tau > 0, same c.o.v.
383 */
384 p[0] = log(lambda);
385 p[1] = log(tau);
386
387 /* pass problem to the optimizer */
388 status = esl_min_ConjugateGradientDescent(NULL, p, 2,
389 &wei_func, NULL,
390 (void *)(&data), &fx, NULL);
391 *ret_mu = mu;
392 *ret_lambda = exp(p[0]);
393 *ret_tau = exp(p[1]);
394 return status;
395 }
396
397
398 /*****************************************************************
399 * 6. ML fitting to binned data
400 *****************************************************************/
401
402
403 struct wei_binned_data {
404 ESL_HISTOGRAM *h; /* contains the binned observed data */
405 double mu; /* mu is considered to be known, not fitted */
406 };
407
408 /* wei_binned_func():
409 * Returns the negative log likelihood of a binned data sample,
410 * in the API of the conjugate gradient descent optimizer in esl_minimizer.
411 */
412 static double
wei_binned_func(double * p,int nparam,void * dptr)413 wei_binned_func(double *p, int nparam, void *dptr)
414 {
415 struct wei_binned_data *data = (struct wei_binned_data *) dptr;
416 ESL_HISTOGRAM *h = data->h;
417 double lambda, tau;
418 double logL;
419 double ai,bi;
420 int i;
421 double tmp;
422
423 /* Unpack what the optimizer gave us.
424 */
425 lambda = exp(p[0]); /* see below for c.o.v. notes */
426 tau = exp(p[1]);
427
428 logL = 0.;
429 for (i = h->cmin; i <= h->imax; i++)
430 {
431 if (h->obs[i] == 0) continue;
432
433 ai = esl_histogram_Bin2LBound(h,i);
434 bi = esl_histogram_Bin2UBound(h,i);
435 if (ai < data->mu) ai = data->mu;
436
437 tmp = esl_wei_cdf(bi, data->mu, lambda, tau) -
438 esl_wei_cdf(ai, data->mu, lambda, tau);
439
440 /* for cdf~1.0, numerical roundoff error can create tmp<0 by a
441 * teensy amount; tolerate that, but catch anything worse */
442 ESL_DASSERT1( (tmp + 1e-7 > 0.));
443 if (tmp <= 0.) return eslINFINITY;
444
445 logL += h->obs[i] * log(tmp);
446 }
447 return -logL; /* goal: minimize NLL */
448 }
449
450 /* Function: esl_wei_FitCompleteBinned()
451 *
452 * Purpose: Given a histogram <g> with binned observations, where each
453 * bin i holds some number of observed samples x with values from
454 * lower bound l to upper bound u (that is, $l < x \leq u$), and
455 * <mu>, the known offset (minimum value) of the distribution;
456 * return maximum likelihood parameters <ret_lambda>
457 * and <ret_tau>.
458 *
459 * Args: x - complete GEV-distributed data [0..n-1]
460 * n - number of samples in <x>
461 * ret_mu - lower bound of the distribution (all x_i > mu)
462 * ret_lambda - RETURN: maximum likelihood estimate of lambda
463 * ret_tau - RETURN: maximum likelihood estimate of tau
464 *
465 * Returns: <eslOK> on success.
466 *
467 * Throws: <eslENOHALT> if the fit doesn't converge.
468 *
469 * Xref: STL9/136-137
470 */
471 int
esl_wei_FitCompleteBinned(ESL_HISTOGRAM * h,double * ret_mu,double * ret_lambda,double * ret_tau)472 esl_wei_FitCompleteBinned(ESL_HISTOGRAM *h, double *ret_mu,
473 double *ret_lambda, double *ret_tau)
474 {
475 struct wei_binned_data data;
476 double p[2]; /* parameter vector */
477 double mean;
478 double mu, lambda, tau; /* initial param guesses */
479 double fx; /* f(x) at minimum; currently unused */
480 int i;
481 double ai;
482 int status;
483
484 /* Set the fixed mu.
485 * Make a good initial guess of lambda, based on exponential fit.
486 * Choose an arbitrary tau.
487 */
488 if (h->is_tailfit) mu = h->phi; /* all x > mu in this case */
489 else if (h->is_rounded) mu = esl_histogram_Bin2LBound(h, h->imin);
490 else mu = h->xmin;
491
492 mean = 0.;
493 for (i = h->cmin; i <= h->imax; i++)
494 {
495 ai = esl_histogram_Bin2LBound(h, i);
496 ai += 0.5*h->w; /* midpoint in bin */
497 mean += (double)h->obs[i] * ai;
498 }
499 mean /= h->No;
500 lambda = 1 / (mean - mu);
501
502 tau = 0.9;
503
504 /* load the data structure */
505 data.h = h;
506 data.mu = mu;
507
508 /* Change of variables;
509 * lambda > 0, so c.o.v. lambda = exp^w, w = log(lambda);
510 * tau > 0, same c.o.v.
511 */
512 p[0] = log(lambda);
513 p[1] = log(tau);
514
515 /* pass problem to the optimizer
516 */
517 status = esl_min_ConjugateGradientDescent(NULL, p, 2,
518 &wei_binned_func, NULL,
519 (void *)(&data), &fx, NULL);
520 *ret_mu = mu;
521 *ret_lambda = exp(p[0]);
522 *ret_tau = exp(p[1]);
523 return status;
524 }
525
526 /*--------------------------- end fitting ----------------------------------*/
527
528
529
530
531 /****************************************************************************
532 * 7. Test driver
533 ****************************************************************************/
534 #ifdef eslWEIBULL_TESTDRIVE
535 /* Compile:
536 gcc -g -Wall -I. -I ~/src/easel -L ~/src/easel -o test -DeslWEIBULL_TESTDRIVE\
537 esl_weibull.c -leasel -lm
538 */
539 #include <stdio.h>
540 #include <stdlib.h>
541 #include <string.h>
542
543 #include "easel.h"
544 #include "esl_getopts.h"
545 #include "esl_random.h"
546 #include "esl_histogram.h"
547 #include "esl_weibull.h"
548
549 static ESL_OPTIONS options[] = {
550 /* name type default env range toggles reqs incomp help docgroup*/
551 { "-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show brief help on version and usage", 0 },
552 { "-l", eslARG_REAL, "1.0", NULL, NULL, NULL, NULL, NULL, "set slope of sampled variates (lambda parameter) to <x> ", 0 },
553 { "-m", eslARG_REAL, "10.0", NULL, NULL, NULL, NULL, NULL, "set location of sampled variates (mu parameter) to <x>", 0 },
554 { "-n", eslARG_INT, "10000", NULL, NULL, NULL, NULL, NULL, "set # of sampled variates to <n>", 0 },
555 { "-o", eslARG_OUTFILE, NULL, NULL, NULL, NULL, NULL, NULL, "output histogram to file <f>", 0 },
556 { "-s", eslARG_INT, "0", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0 },
557 { "-t", eslARG_REAL, "0.7", NULL, NULL, NULL, NULL, NULL, "set shape of sampled variates (tau parameter) to <x>", 0 },
558 { "-v", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "be more verbose in output", 0 },
559 { "-w", eslARG_REAL, "0.1", NULL, NULL, NULL, NULL, NULL, "set width of histogram bins to <x>", 0 },
560 { "--cdf", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of cumulative distribution", 0 },
561 { "--logcdf", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of log cumulative distribution", 0 },
562 { "--pdf", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of probability density", 0 },
563 { "--logpdf", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of log probability density", 0 },
564 { "--surv", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of survival P(s>x)", 0 },
565 { "--logsurv", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "dump plot of log survival, log(P(s>x))", 0 },
566 { "--xL", eslARG_REAL, NULL, NULL, NULL, NULL, NULL, NULL, "set minimum x-axis value on dumped plots to <x>", 0 },
567 { "--xH", eslARG_REAL, NULL, NULL, NULL, NULL, NULL, NULL, "set maximum x-axis value on dumped plots to <x>", 0 },
568 { "--xS", eslARG_REAL, NULL, NULL, NULL, NULL, NULL, NULL, "set x-axis increment value on dumped plots to <x>", 0 },
569 { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
570 };
571 static char usage[] = "[-options]";
572 static char banner[] = "test driver for Easel's Weibull distribution module";
573
574 int
main(int argc,char ** argv)575 main(int argc, char **argv)
576 {
577 ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
578 ESL_RANDOMNESS *rng = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
579 double mu = esl_opt_GetReal (go, "-m");
580 double lambda = esl_opt_GetReal (go, "-l");
581 double tau = esl_opt_GetReal (go, "-t");
582 int n = esl_opt_GetInteger(go, "-n");
583 double binwidth = esl_opt_GetReal (go, "-w");
584 int plot_cdf = esl_opt_GetBoolean(go, "--cdf");
585 int plot_logcdf = esl_opt_GetBoolean(go, "--logcdf");
586 int plot_pdf = esl_opt_GetBoolean(go, "--pdf");
587 int plot_logpdf = esl_opt_GetBoolean(go, "--logpdf");
588 int plot_surv = esl_opt_GetBoolean(go, "--surv");
589 int plot_logsurv = esl_opt_GetBoolean(go, "--logsurv");
590 int be_verbose = esl_opt_GetBoolean(go, "-v");
591 char *plotfile = esl_opt_GetString (go, "-o");
592 ESL_HISTOGRAM *h = NULL;
593 int xmin_set = esl_opt_IsOn(go, "--xL");
594 double xmin = xmin_set ? esl_opt_GetReal(go, "--xL") : mu;
595 int xmax_set = esl_opt_IsOn(go, "--xH");
596 double xmax = xmax_set ? esl_opt_GetReal(go, "--xH") : mu+40*(1./lambda);
597 int xstep_set = esl_opt_IsOn(go, "--xH");
598 double xstep = xstep_set ? esl_opt_GetReal(go, "--xS") : 0.1;
599 FILE *pfp = stdout;
600 double emu, elambda, etau;
601 int i;
602 double x;
603 double *data;
604 int ndata;
605
606 fprintf(stderr, "## %s\n", argv[0]);
607 fprintf(stderr, "# rng seed = %" PRIu32 "\n", esl_randomness_GetSeed(rng));
608
609 if (be_verbose) printf("Parametric: mu = %f lambda = %f tau = %f\n", mu, lambda, tau);
610
611 h = esl_histogram_CreateFull(mu, 100., binwidth);
612 if (plotfile && (pfp = fopen(plotfile, "w")) == NULL) ESL_EXCEPTION(eslFAIL, "Failed to open plotfile");
613
614 for (i = 0; i < n; i++)
615 {
616 x = esl_wei_Sample(rng, mu, lambda, tau);
617 esl_histogram_Add(h, x);
618 }
619 esl_histogram_GetData(h, &data, &ndata);
620
621 esl_wei_FitComplete(data, ndata, &emu, &elambda, &etau);
622 if (be_verbose) printf("Complete data fit: mu = %f lambda = %f tau = %f\n", emu, elambda, etau);
623 if (fabs( (emu-mu)/mu ) > 0.01) ESL_EXCEPTION(eslFAIL, "Error in (complete) fitted mu > 1%\n");
624 if (fabs( (elambda-lambda)/lambda ) > 0.10) ESL_EXCEPTION(eslFAIL, "Error in (complete) fitted lambda > 10%\n");
625 if (fabs( (etau-tau)/tau ) > 0.10) ESL_EXCEPTION(eslFAIL, "Error in (complete) fitted tau > 10%\n");
626
627 esl_wei_FitCompleteBinned(h, &emu, &elambda, &etau);
628 if (be_verbose) printf("Binned data fit: mu = %f lambda = %f tau = %f\n", emu, elambda, etau);
629 if (fabs( (emu-mu)/mu ) > 0.01) ESL_EXCEPTION(eslFAIL, "Error in (binned) fitted mu > 1%\n");
630 if (fabs( (elambda-lambda)/lambda ) > 0.10) ESL_EXCEPTION(eslFAIL, "Error in (binned) fitted lambda > 10%\n");
631 if (fabs( (etau-tau)/tau ) > 0.10) ESL_EXCEPTION(eslFAIL, "Error in (binned) fitted lambda > 10%\n");
632
633 if (plot_pdf) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_pdf, xmin, xmax, xstep);
634 if (plot_logpdf) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_logpdf, xmin, xmax, xstep);
635 if (plot_cdf) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_cdf, xmin, xmax, xstep);
636 if (plot_logcdf) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_logcdf, xmin, xmax, xstep);
637 if (plot_surv) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_surv, xmin, xmax, xstep);
638 if (plot_logsurv) esl_wei_Plot(pfp, mu, lambda, tau, &esl_wei_logsurv, xmin, xmax, xstep);
639
640 if (plotfile) fclose(pfp);
641 esl_histogram_Destroy(h);
642 esl_randomness_Destroy(rng);
643 esl_getopts_Destroy(go);
644
645 fprintf(stderr, "# status = ok\n");
646 return 0;
647 }
648 #endif /*eslWEIBULL_TESTDRIVE*/
649
650 /****************************************************************************
651 * 8. Example
652 ****************************************************************************/
653 #ifdef eslWEIBULL_EXAMPLE
654 /*::cexcerpt::wei_example::begin::*/
655 #include <stdio.h>
656 #include "easel.h"
657 #include "esl_random.h"
658 #include "esl_histogram.h"
659 #include "esl_weibull.h"
660
661 int
main(int argc,char ** argv)662 main(int argc, char **argv)
663 {
664 double mu = -2.1;
665 double lambda = 1.0;
666 double tau = 0.8;
667 ESL_HISTOGRAM *h = esl_histogram_CreateFull(mu, 100., 0.1);
668 ESL_RANDOMNESS *r = esl_randomness_Create(0);
669 int n = 10000;
670 double emu, elambda, etau;
671 double *data;
672 int ndata;
673 double x;
674 int i;
675
676 for (i = 0; i < n; i++)
677 {
678 x = esl_wei_Sample(r, mu, lambda, tau);
679 esl_histogram_Add(h, x);
680 }
681 esl_histogram_GetData(h, &data, &ndata);
682
683 /* Plot the empirical (sampled) and expected survivals */
684 esl_histogram_PlotSurvival(stdout, h);
685 esl_wei_Plot(stdout, mu, lambda, tau,
686 &esl_wei_surv, h->xmin, h->xmax, 0.1);
687
688 /* ML fit to complete data, and plot fitted survival curve */
689 esl_wei_FitComplete(data, ndata, &emu, &elambda, &etau);
690 esl_wei_Plot(stdout, emu, elambda, etau,
691 &esl_wei_surv, h->xmin, h->xmax, 0.1);
692
693 /* ML fit to binned data, plot fitted survival curve */
694 esl_wei_FitCompleteBinned(h, &emu, &elambda, &etau);
695 esl_wei_Plot(stdout, emu, elambda, etau,
696 &esl_wei_surv, h->xmin, h->xmax, 0.1);
697
698 esl_randomness_Destroy(r);
699 esl_histogram_Destroy(h);
700 return 0;
701 }
702 /*::cexcerpt::wei_example::end::*/
703 #endif /*eslWEIBULL_EXAMPLE*/
704
705
706
707