1 /* Statistical routines for stretched exponential distributions.
2 *
3 * Contents:
4 * 1. Evaluating densities and distributions
5 * 2. Generic API routines: for general interface w/ histogram module
6 * 3. Dumping plots for 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 * Xrefs:
14 * STL9/146 : original implementation
15 *
16 * To-do:
17 * - Fit*() functions should return eslEINVAL on n=0, eslENORESULT
18 * on failure due to small n. Compare esl_gumbel. xref J12/93.
19 * SRE, Wed Nov 27 11:07:44 2013
20 */
21 #include "esl_config.h"
22
23 #include <stdio.h>
24 #include <math.h>
25
26 #include "easel.h"
27 #include "esl_histogram.h"
28 #include "esl_minimizer.h"
29 #include "esl_random.h"
30 #include "esl_stats.h"
31 #include "esl_vectorops.h"
32 #include "esl_stretchexp.h"
33
34
35 /****************************************************************************
36 * 1. Evaluating densities and distributions
37 ****************************************************************************/
38 /* mu <= x < infinity
39 * [x=mu is no problem, but watch out for evaluating log(0) when it is]
40 * lambda > 0
41 * tau > 0 [fat tailed when tau < 1; thin when tau > 1; exponential when tau = 1]
42 */
43
44 /* Function: esl_sxp_pdf()
45 *
46 * Purpose: Calculates the probability density function for the
47 * stretched exponential pdf, $P(X=x)$, given
48 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
49 */
50 double
esl_sxp_pdf(double x,double mu,double lambda,double tau)51 esl_sxp_pdf(double x, double mu, double lambda, double tau)
52 {
53 double y = lambda * (x-mu);
54 double val;
55 double gt;
56
57 if (x < mu) return 0.;
58 esl_stats_LogGamma(1/tau, >);
59
60 if (x == mu) val = (lambda * tau / exp(gt));
61 else val = (lambda * tau / exp(gt)) * exp(- exp(tau * log(y)));
62
63 return val;
64 }
65
66 /* Function: esl_sxp_logpdf()
67 *
68 * Purpose: Calculates the log probability density function for the
69 * stretched exponential pdf, $\log P(X=x)$, given
70 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
71 */
72 double
esl_sxp_logpdf(double x,double mu,double lambda,double tau)73 esl_sxp_logpdf(double x, double mu, double lambda, double tau)
74 {
75 double y = lambda * (x-mu);
76 double gt;
77 double val;
78
79 if (x < mu) return -eslINFINITY;
80 esl_stats_LogGamma(1/tau, >);
81
82 if (x == mu) val = log(lambda) + log(tau) - gt;
83 else val = log(lambda) + log(tau) - gt - exp(tau*log(y));
84 return val;
85 }
86
87 /* Function: esl_sxp_cdf()
88 *
89 * Purpose: Calculates the cumulative distribution function for the
90 * stretched exponential pdf, $P(X \leq x)$, given
91 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
92 */
93 double
esl_sxp_cdf(double x,double mu,double lambda,double tau)94 esl_sxp_cdf(double x, double mu, double lambda, double tau)
95 {
96 double y = lambda * (x-mu);
97 double val;
98
99 if (x <= mu) return 0.;
100 esl_stats_IncompleteGamma(1/tau, exp(tau * log(y)), &val, NULL);
101
102 ESL_DASSERT1 (( !isnan(val)));
103 return val;
104 }
105
106 /* Function: esl_sxp_logcdf()
107 *
108 * Purpose: Calculates the log of the cumulative distribution function for the
109 * stretched exponential pdf, $\log P(X \leq x)$, given
110 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
111 */
112 double
esl_sxp_logcdf(double x,double mu,double lambda,double tau)113 esl_sxp_logcdf(double x, double mu, double lambda, double tau)
114 {
115 double y = lambda * (x-mu);
116 double val;
117
118 if (x <= mu) return -eslINFINITY;
119 esl_stats_IncompleteGamma(1./tau, exp(tau * log(y)), &val, NULL);
120 return log(val);
121 }
122
123 /* Function: esl_sxp_surv()
124 *
125 * Purpose: Calculates the survival function for the
126 * stretched exponential pdf, $P(X > x)$, given
127 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
128 */
129 double
esl_sxp_surv(double x,double mu,double lambda,double tau)130 esl_sxp_surv(double x, double mu, double lambda, double tau)
131 {
132 double y = lambda * (x-mu);
133 double val;
134
135 if (x <= mu) return 1.0;
136
137 esl_stats_IncompleteGamma(1./tau, exp(tau * log(y)), NULL, &val);
138 return val;
139 }
140
141 /* Function: esl_sxp_logsurv()
142 *
143 * Purpose: Calculates the log survival function for the
144 * stretched exponential pdf, $\log P(X > x)$, given
145 * quantile <x>, offset <mu>, and parameters <lambda> and <tau>.
146 */
147 double
esl_sxp_logsurv(double x,double mu,double lambda,double tau)148 esl_sxp_logsurv(double x, double mu, double lambda, double tau)
149 {
150 double y = lambda * (x-mu);
151 double val;
152
153 if (x <= mu) return 0.0;
154
155 esl_stats_IncompleteGamma(1./tau, exp(tau * log(y)), NULL, &val);
156 return log(val);
157 }
158
159 /* Function: esl_sxp_invcdf()
160 *
161 * Purpose: Calculates the inverse CDF for a stretched exponential
162 * with parameters <mu>, <lambda>, and <tau>, returning
163 * the quantile <x> at which the CDF is <p>.
164 *
165 * The inverse CDF of the stretched exponential has no
166 * analytical expression as far as I'm aware. The calculation
167 * here is a computationally expensive, brute force bisection
168 * search in <x> using the CDF function. It will suffice for
169 * a small number of calls (for plotting applications, for example),
170 * but it is not sufficient for a large number of calls.
171 */
172 double
esl_sxp_invcdf(double p,double mu,double lambda,double tau)173 esl_sxp_invcdf(double p, double mu, double lambda, double tau)
174 {
175 double x1, x2, xm; /* low, high guesses at x */
176 double f2, fm;
177 double tol = 1e-6;
178
179 x1 = mu;
180 x2 = mu + 1.;
181 do { /* bracket */
182 x2 = x2 + 2.*(x2-x1);
183 f2 = esl_sxp_cdf(x2, mu, lambda, tau);
184 } while (f2 < p);
185
186 do { /* bisection */
187 xm = (x1+x2) / 2.;
188 fm = esl_sxp_cdf(xm, mu, lambda, tau);
189
190 if (fm > p) x2 = xm;
191 else if (fm < p) x1 = xm;
192 else return xm; /* unlikely case of fm==cdf */
193 } while ( (x2-x1)/(x1+x2-2*mu) > tol);
194
195 xm = (x1+x2) / 2.;
196 return xm;
197 }
198 /*-------------------- end densities & distributions ------------------------*/
199
200
201
202
203 /****************************************************************************
204 * 2. Generic API routines: for general interface w/ histogram module
205 ****************************************************************************/
206
207 /* Function: esl_sxp_generic_pdf()
208 *
209 * Purpose: Generic-API wrapper around <esl_sxp_pdf()>, taking
210 * a void ptr to a double array containing $\mu$, $\lambda$,
211 * $\tau$ parameters.
212 */
213 double
esl_sxp_generic_pdf(double x,void * params)214 esl_sxp_generic_pdf(double x, void *params)
215 {
216 double *p = (double *) params;
217 return esl_sxp_pdf(x, p[0], p[1], p[2]);
218 }
219
220 /* Function: esl_sxp_generic_cdf()
221 *
222 * Purpose: Generic-API wrapper around <esl_sxp_cdf()>, taking
223 * a void ptr to a double array containing $\mu$, $\lambda$,
224 * $\tau$ parameters.
225 */
226 double
esl_sxp_generic_cdf(double x,void * params)227 esl_sxp_generic_cdf(double x, void *params)
228 {
229 double *p = (double *) params;
230 return esl_sxp_cdf(x, p[0], p[1], p[2]);
231 }
232
233 /* Function: esl_sxp_generic_surv()
234 *
235 * Purpose: Generic-API wrapper around <esl_sxp_surv()>, taking
236 * a void ptr to a double array containing $\mu$, $\lambda$,
237 * $\tau$ parameters.
238 */
239 double
esl_sxp_generic_surv(double x,void * params)240 esl_sxp_generic_surv(double x, void *params)
241 {
242 double *p = (double *) params;
243 return esl_sxp_surv(x, p[0], p[1], p[2]);
244 }
245
246 /* Function: esl_sxp_generic_invcdf()
247 *
248 * Purpose: Generic-API wrapper around <esl_sxp_invcdf()>, taking
249 * a void ptr to a double array containing $\mu$, $\lambda$,
250 * $\tau$ parameters.
251 */
252 double
esl_sxp_generic_invcdf(double p,void * params)253 esl_sxp_generic_invcdf(double p, void *params)
254 {
255 double *v = (double *) params;
256 return esl_sxp_invcdf(p, v[0], v[1], v[2]);
257 }
258 /*------------------------ end generic API ---------------------------------*/
259
260
261
262 /****************************************************************************
263 * 3. Dumping plots for files
264 ****************************************************************************/
265
266 /* Function: esl_sxp_Plot()
267 *
268 * Purpose: Plot some stretched exponential function <func> (for instance,
269 * <esl_sxp_pdf()>) for parameters <mu>, <lambda>, and <tau>, for
270 * a range of quantiles x from <xmin> to <xmax> in steps of <xstep>;
271 * output to an open stream <fp> in xmgrace XY input format.
272 *
273 * Returns: <eslOK> on success.
274 *
275 * Throws: <eslEWRITE> on any system write error, such as filled disk.
276 */
277 int
esl_sxp_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)278 esl_sxp_Plot(FILE *fp, double mu, double lambda, double tau,
279 double (*func)(double x, double mu, double lambda, double tau),
280 double xmin, double xmax, double xstep)
281 {
282 double x;
283 for (x = xmin; x <= xmax; x += xstep)
284 if (fprintf(fp, "%f\t%g\n", x, (*func)(x, mu, lambda, tau)) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "stretchexp plot write failed");
285 if (fprintf(fp, "&\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "stretchexp plot write failed");
286 return eslOK;
287 }
288 /*-------------------- end plot dumping routines ---------------------------*/
289
290
291
292
293 /****************************************************************************
294 * 4. Sampling
295 ****************************************************************************/
296
297 /* Function: esl_sxp_Sample()
298 *
299 * Purpose: Sample a stretched exponential random variate,
300 * by a change of variable from a Gamma sample.
301 */
302 double
esl_sxp_Sample(ESL_RANDOMNESS * r,double mu,double lambda,double tau)303 esl_sxp_Sample(ESL_RANDOMNESS *r, double mu, double lambda, double tau)
304 {
305 double t,x;
306
307 t = esl_rnd_Gamma(r, 1./tau);
308 x = mu + 1./lambda * exp(1./tau * log(t));
309 return x;
310 }
311 /*--------------------------- end sampling ---------------------------------*/
312
313
314
315 /****************************************************************************
316 * 5. ML fitting to complete data
317 ****************************************************************************/
318
319 /* This structure is used to sneak the data into minimizer's generic
320 * (void *) API for all aux data
321 */
322 struct sxp_data {
323 double *x;
324 int n;
325 double mu;
326 };
327
328 static double
sxp_complete_func(double * p,int np,void * dptr)329 sxp_complete_func(double *p, int np, void *dptr)
330 {
331 struct sxp_data *data = (struct sxp_data *) dptr;
332 double lambda, tau;
333 double logL = 0.;
334 int i;
335
336 lambda = exp(p[0]);
337 tau = exp(p[1]);
338
339 for (i = 0; i < data->n; i++)
340 logL += esl_sxp_logpdf(data->x[i], data->mu, lambda, tau);
341 return -logL;
342 }
343
344 /* Function: esl_sxp_FitComplete()
345 *
346 * Purpose: Given a vector of <n> observed data samples <x[]>,
347 * find maximum likelihood parameters by conjugate gradient
348 * descent optimization.
349 */
350 int
esl_sxp_FitComplete(double * x,int n,double * ret_mu,double * ret_lambda,double * ret_tau)351 esl_sxp_FitComplete(double *x, int n,
352 double *ret_mu, double *ret_lambda, double *ret_tau)
353
354 {
355 struct sxp_data data;
356 double p[2];
357 double mu, tau, lambda;
358 double mean;
359 double fx;
360 int status;
361
362 /* initial guesses; mu is definitely = minimum x,
363 * and just use arbitrary #'s to init lambda, tau
364 */
365 mu = esl_vec_DMin(x, n);
366 esl_stats_DMean(x, n, &mean, NULL);
367 lambda = 1 / (mean - mu);
368 tau = 0.9;
369
370
371 /* load data structure, param vector, and step vector */
372 data.x = x;
373 data.n = n;
374 data.mu = mu;
375 p[0] = log(lambda);
376 p[1] = log(tau);
377
378 /* hand it off */
379 status = esl_min_ConjugateGradientDescent(NULL, p, 2,
380 &sxp_complete_func, NULL,
381 (void *) (&data), &fx, NULL);
382 *ret_mu = mu;
383 *ret_lambda = exp(p[0]);
384 *ret_tau = exp(p[1]);
385 return status;
386 }
387
388
389 /****************************************************************************
390 * 6. ML fitting to binned data
391 ****************************************************************************/
392
393 struct sxp_binned_data {
394 ESL_HISTOGRAM *g; /* contains the binned data */
395 double mu; /* mu is not a learnable param */
396 };
397
398 static double
sxp_complete_binned_func(double * p,int np,void * dptr)399 sxp_complete_binned_func(double *p, int np, void *dptr)
400 {
401 struct sxp_binned_data *data = (struct sxp_binned_data *) dptr;
402 ESL_HISTOGRAM *g = data->g;
403 double logL = 0.;
404 double ai, bi; /* lower, upper bounds on bin */
405 double lambda, tau;
406 int i;
407 double tmp;
408
409 lambda = exp(p[0]);
410 tau = exp(p[1]);
411
412 ESL_DASSERT1(( ! isnan(lambda) ));
413 ESL_DASSERT1(( ! isnan(tau) ));
414
415 for (i = g->cmin; i <= g->imax; i++) /* for each occupied bin */
416 {
417 if (g->obs[i] == 0) continue;
418
419 ai = esl_histogram_Bin2LBound(g, i);
420 bi = esl_histogram_Bin2UBound(g, i);
421 if (ai < data->mu) ai = data->mu; /* careful at leftmost bound */
422
423 tmp = esl_sxp_cdf(bi, data->mu, lambda, tau) -
424 esl_sxp_cdf(ai, data->mu, lambda, tau);
425 if (tmp == 0.) return eslINFINITY;
426 logL += g->obs[i] * log(tmp);
427 }
428 return -logL; /* minimizing NLL */
429 }
430
431 /* Function: esl_sxp_FitCompleteBinned()
432 *
433 * Purpose: Given a histogram <g> with binned observations, where each
434 * bin i holds some number of observed samples x with values from
435 * lower bound l to upper bound u (that is, $l < x \leq u$);
436 * find maximum likelihood parameters mu, lambda, tau by conjugate
437 * gradient descent optimization.
438 */
439 int
esl_sxp_FitCompleteBinned(ESL_HISTOGRAM * g,double * ret_mu,double * ret_lambda,double * ret_tau)440 esl_sxp_FitCompleteBinned(ESL_HISTOGRAM *g,
441 double *ret_mu, double *ret_lambda, double *ret_tau)
442
443 {
444 struct sxp_binned_data data;
445 double p[2];
446 double mu, tau, lambda;
447 double fx;
448 double ai, mean;
449 int i;
450 int status;
451
452 /* Set the fixed mu.
453 * Make a good initial guess of lambda, based on exponential fit.
454 * Choose an arbitrary tau.
455 */
456 if (g->is_tailfit) mu = g->phi; /* all x > mu in this case */
457 else if (g->is_rounded) mu = esl_histogram_Bin2LBound(g, g->imin);
458 else mu = g->xmin;
459
460 mean = 0.;
461 for (i = g->cmin; i <= g->imax; i++)
462 {
463 ai = esl_histogram_Bin2LBound(g, i);
464 ai += 0.5*g->w; /* midpoint in bin */
465 mean += (double)g->obs[i] * ai;
466 }
467 mean /= g->No;
468 lambda = 1 / (mean - mu);
469 tau = 0.9;
470
471 /* load data structure, param vector, and step vector */
472 data.g = g;
473 data.mu = mu;
474 p[0] = log(lambda);
475 p[1] = log(tau);
476
477 /* hand it off */
478 status = esl_min_ConjugateGradientDescent(NULL, p, 2,
479 &sxp_complete_binned_func, NULL,
480 (void *) (&data), &fx, NULL);
481 *ret_mu = mu;
482 *ret_lambda = exp(p[0]);
483 *ret_tau = exp(p[1]);
484 return status;
485 }
486
487
488
489
490 /****************************************************************************
491 * 7. Test driver
492 ****************************************************************************/
493 #ifdef eslSTRETCHEXP_TESTDRIVE
494 /* gcc -g -Wall -I. -L . -o stretchexp_utest -DeslSTRETCHEXP_TESTDRIVE esl_stretchexp.c -leasel -lm
495 */
496 #include "esl_config.h"
497
498 #include <stdio.h>
499 #include <stdlib.h>
500 #include <string.h>
501
502 #include "easel.h"
503 #include "esl_getopts.h"
504 #include "esl_random.h"
505 #include "esl_histogram.h"
506 #include "esl_stretchexp.h"
507
508 static ESL_OPTIONS options[] = {
509 /* name type default env range togs reqs incomp help docgrp */
510 {"-h", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show help and usage", 0},
511 {"-l", eslARG_REAL, "1.0", NULL,"x>0", NULL, NULL, NULL, "set lambda param to <x>", 0},
512 {"-m", eslARG_REAL, "10.0", NULL, NULL, NULL, NULL, NULL, "set mu param to <x>", 0},
513 {"-n", eslARG_INT, "10000", NULL,"n>0", NULL, NULL, NULL, "set number of samples to <n>", 0},
514 {"-o", eslARG_OUTFILE, NULL, NULL, NULL, NULL, NULL, NULL, "save plots to file <f>", 0},
515 {"-s", eslARG_INT, "42", NULL, NULL, NULL, NULL, NULL, "set random number seed to <n>", 0},
516 {"-t", eslARG_REAL, "0.7", NULL,"x>0", NULL, NULL, NULL, "set tau param to <x>", 0},
517 {"-v", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "show verbose output", 0},
518 {"-w", eslARG_REAL, "0.1", NULL,"x>0", NULL, NULL, NULL, "set width of histogram bins to <x>",0},
519 {"--C", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot CDF", 0},
520 {"--LC",eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot log CDF", 0},
521 {"--P", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot PDF", 0},
522 {"--LP",eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot log PDF", 0},
523 {"--S", eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot survival", 0},
524 {"--LS",eslARG_NONE, FALSE, NULL, NULL, NULL, NULL, NULL, "plot log survival", 0},
525 {"--XL",eslARG_NONE, NULL, NULL, NULL, NULL, NULL, NULL, "set xmin for plot axis", 0},
526 {"--XH",eslARG_NONE, NULL, NULL, NULL, NULL, NULL, NULL, "set xmax for plot axis", 0},
527 {"--XS",eslARG_NONE, NULL, NULL, NULL, NULL, NULL, NULL, "set xstep for plot axis", 0},
528 { 0,0,0,0,0,0,0,0,0,0},
529 };
530 static char usage[] = "[-options]";
531 static char banner[] = "test driver for random module";
532
533 int
main(int argc,char ** argv)534 main(int argc, char **argv)
535 {
536 ESL_GETOPTS *go = esl_getopts_CreateDefaultApp(options, 0, argc, argv, banner, usage);
537 double mu = esl_opt_GetReal(go, "-m");
538 double lambda = esl_opt_GetReal(go, "-l");
539 double tau = esl_opt_GetReal(go, "-t");
540 ESL_RANDOMNESS *r = esl_randomness_Create(esl_opt_GetInteger(go, "-s"));
541 ESL_HISTOGRAM *h = esl_histogram_CreateFull(mu, 100., esl_opt_GetReal(go, "-w"));;
542 int n = esl_opt_GetInteger(go, "-n");
543 int be_verbose = esl_opt_GetBoolean(go, "-v");
544 char *plotfile = esl_opt_GetString(go, "-o");
545 FILE *pfp = stdout;
546 int plot_pdf = esl_opt_GetBoolean(go, "--P");
547 int plot_logpdf = esl_opt_GetBoolean(go, "--LP");
548 int plot_cdf = esl_opt_GetBoolean(go, "--C");
549 int plot_logcdf = esl_opt_GetBoolean(go, "--LC");
550 int plot_surv = esl_opt_GetBoolean(go, "--S");
551 int plot_logsurv = esl_opt_GetBoolean(go, "--LS");
552 double xmin = esl_opt_IsOn(go, "--XL") ? esl_opt_GetReal(go, "--XL") : mu;
553 double xmax = esl_opt_IsOn(go, "--XH") ? esl_opt_GetReal(go, "--XH") : mu+40*(1./lambda);
554 double xstep = esl_opt_IsOn(go, "--XS") ? esl_opt_GetReal(go, "--XS") : 0.1;
555 double emu, elambda, etau;
556 int i;
557 double x;
558 double *data;
559 int ndata;
560
561 if (be_verbose)
562 printf("Parametric: mu = %f lambda = %f tau = %f\n", mu, lambda, tau);
563
564 if (plotfile != NULL) {
565 if ((pfp = fopen(plotfile, "w")) == NULL)
566 esl_fatal("Failed to open plotfile");
567 }
568
569 for (i = 0; i < n; i++)
570 {
571 x = esl_sxp_Sample(r, mu, lambda, tau);
572 esl_histogram_Add(h, x);
573 }
574 esl_histogram_GetData(h, &data, &ndata);
575
576 esl_sxp_FitComplete(data, ndata, &emu, &elambda, &etau);
577 if (be_verbose)
578 printf("Complete data fit: mu = %f lambda = %f tau = %f\n",
579 emu, elambda, etau);
580 if (fabs( (emu-mu)/mu ) > 0.01) esl_fatal("Error in (complete) fitted mu > 1%\n");
581 if (fabs( (elambda-lambda)/lambda ) > 0.10) esl_fatal("Error in (complete) fitted lambda > 10%\n");
582 if (fabs( (etau-tau)/tau ) > 0.10) esl_fatal("Error in (complete) fitted tau > 10%\n");
583
584 esl_sxp_FitCompleteBinned(h, &emu, &elambda, &etau);
585 if (be_verbose)
586 printf("Binned data fit: mu = %f lambda = %f tau = %f\n",
587 emu, elambda, etau);
588 if (fabs( (emu-mu)/mu ) > 0.01) esl_fatal("Error in (binned) fitted mu > 1%\n");
589 if (fabs( (elambda-lambda)/lambda ) > 0.10) esl_fatal("Error in (binned) fitted lambda > 10%\n");
590 if (fabs( (etau-tau)/tau ) > 0.10) esl_fatal("Error in (binned) fitted tau > 10%\n");
591
592 if (plot_pdf) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_pdf, xmin, xmax, xstep);
593 if (plot_logpdf) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_logpdf, xmin, xmax, xstep);
594 if (plot_cdf) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_cdf, xmin, xmax, xstep);
595 if (plot_logcdf) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_logcdf, xmin, xmax, xstep);
596 if (plot_surv) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_surv, xmin, xmax, xstep);
597 if (plot_logsurv) esl_sxp_Plot(pfp, mu, lambda, tau, &esl_sxp_logsurv, xmin, xmax, xstep);
598
599 if (plotfile != NULL) fclose(pfp);
600 esl_histogram_Destroy(h);
601 esl_randomness_Destroy(r);
602 esl_getopts_Destroy(go);
603 return 0;
604 }
605 #endif /*eslSTRETCHEXP_TESTDRIVE*/
606
607 /****************************************************************************
608 * Example main()
609 ****************************************************************************/
610 #ifdef eslSTRETCHEXP_EXAMPLE
611 /*::cexcerpt::sxp_example::begin::*/
612 #include <stdio.h>
613 #include "easel.h"
614 #include "esl_random.h"
615 #include "esl_histogram.h"
616 #include "esl_stretchexp.h"
617
618 int
main(int argc,char ** argv)619 main(int argc, char **argv)
620 {
621 double mu = -50.0;
622 double lambda = 2.5;
623 double tau = 0.7;
624 ESL_HISTOGRAM *h = esl_histogram_CreateFull(mu, 100., 0.1);
625 ESL_RANDOMNESS *r = esl_randomness_Create(0);
626 int n = 10000;
627 double *data;
628 int ndata;
629 double emu, elam, etau;
630 int i;
631 double x;
632
633 for (i = 0; i < n; i++)
634 {
635 x = esl_sxp_Sample(r, mu, lambda, tau);
636 esl_histogram_Add(h, x);
637 }
638 esl_histogram_GetData(h, &data, &ndata);
639
640 /* Plot the empirical (sampled) and expected survivals */
641 esl_histogram_PlotSurvival(stdout, h);
642 esl_sxp_Plot(stdout, mu, lambda, tau,
643 &esl_sxp_surv, h->xmin, h->xmax, 0.1);
644
645 /* ML fit to complete data, and plot fitted survival curve */
646 esl_sxp_FitComplete(data, ndata, &emu, &elam, &etau);
647 esl_sxp_Plot(stdout, emu, elam, etau,
648 &esl_sxp_surv, h->xmin, h->xmax, 0.1);
649
650 /* ML fit to binned data, plot fitted survival curve */
651 esl_sxp_FitCompleteBinned(h, &emu, &elam, &etau);
652 esl_sxp_Plot(stdout, emu, elam, etau,
653 &esl_sxp_surv, h->xmin, h->xmax, 0.1);
654
655 esl_randomness_Destroy(r);
656 esl_histogram_Destroy(h);
657 return 0;
658 }
659 /*::cexcerpt::sxp_example::end::*/
660 #endif /*eslSTRETCHEXP_EXAMPLE*/
661
662