1 /* Statistical routines for exponential 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. Routines for dumping plots for files
7 * 4. Routines for sampling
8 * 5. Maximum likelihood fitting
9 * 6. Stats driver
10 * 7. Unit tests
11 * 8. Test driver
12 * 9. Example
13 *
14 * xref STL9/138
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:03:07 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_random.h"
29 #include "esl_stats.h"
30
31 #include "esl_exponential.h"
32
33
34 /****************************************************************************
35 * 1. Routines for evaluating densities and distributions
36 ****************************************************************************/
37 /* lambda > 0
38 * mu <= x < infinity
39 *
40 * watch out:
41 * - any lambda > 0 is valid... including infinity. Fitting code
42 * may try to test such lambdas, and it must get back valid numbers,
43 * never an NaN, or it will fail. IEEE754 allows us
44 * to calculate log(inf) = inf, exp(-inf) = 0, and exp(inf) = inf.
45 * But inf-inf = NaN, so Don't Do That.
46 */
47
48 /* Function: esl_exp_pdf()
49 *
50 * Purpose: Calculates the probability density function for the
51 * exponential, $P(X=x)$, given value <x>, offset <mu>,
52 * and decay parameter <lambda>.
53 */
54 double
esl_exp_pdf(double x,double mu,double lambda)55 esl_exp_pdf(double x, double mu, double lambda)
56 {
57 if (x < mu) return 0.;
58 return (lambda * exp(-lambda*(x-mu)));
59 }
60
61 /* Function: esl_exp_logpdf()
62 *
63 * Purpose: Calculates the log probability density function for the
64 * exponential, $P(X=x)$, given value <x>, offset <mu>,
65 * and decay parameter <lambda>.
66 */
67 double
esl_exp_logpdf(double x,double mu,double lambda)68 esl_exp_logpdf(double x, double mu, double lambda)
69 {
70 if (x < mu) return -eslINFINITY;
71
72 if (lambda == eslINFINITY)
73 { /* limit as lambda->inf: avoid inf-inf! */
74 if (x == mu) return eslINFINITY;
75 else return -eslINFINITY;
76 }
77 return (log(lambda) - lambda*(x-mu));
78 }
79
80 /* Function: esl_exp_cdf()
81 *
82 * Purpose: Calculates the cumulative distribution function for the
83 * exponential, $P(X \leq x)$, given value <x>, offset <mu>,
84 * and decay parameter <lambda>.
85 */
86 double
esl_exp_cdf(double x,double mu,double lambda)87 esl_exp_cdf(double x, double mu, double lambda)
88 {
89 double y = lambda*(x-mu); /* y>=0 because lambda>0 and x>=mu */
90
91 if (x < mu) return 0.;
92
93 /* 1-e^-y ~ y for small |y| */
94 if (y < eslSMALLX1) return y;
95 else return 1 - exp(-y);
96 }
97
98 /* Function: esl_exp_logcdf()
99 *
100 * Purpose: Calculates the log of the cumulative distribution function
101 * for the exponential, $log P(X \leq x)$, given value <x>,
102 * offset <mu>, and decay parameter <lambda>.
103 */
104 double
esl_exp_logcdf(double x,double mu,double lambda)105 esl_exp_logcdf(double x, double mu, double lambda)
106 {
107 double y = lambda * (x-mu);
108 double ey = exp(-y);
109
110 if (x < mu) return -eslINFINITY;
111
112 /* When y is small, 1-e^-y = y, so answer is log(y);
113 * when y is large, exp(-y) is small, log(1-exp(-y)) = -exp(-y).
114 */
115 if (y == 0) return -eslINFINITY; /* don't allow NaN */
116 else if (y < eslSMALLX1) return log(y);
117 else if (ey < eslSMALLX1) return -ey;
118 else return log(1-ey);
119 }
120
121 /* Function: esl_exp_surv()
122 *
123 * Purpose: Calculates the survivor function, $P(X>x)$ (that is, 1-CDF,
124 * the right tail probability mass) for an exponential distribution,
125 * given value <x>, offset <mu>, and decay parameter <lambda>.
126 */
127 double
esl_exp_surv(double x,double mu,double lambda)128 esl_exp_surv(double x, double mu, double lambda)
129 {
130 if (x < mu) return 1.0;
131 return exp(-lambda * (x-mu));
132 }
133
134 /* Function: esl_exp_logsurv()
135 *
136 * Purpose: Calculates the log survivor function, $\log P(X>x)$ (that is,
137 * log(1-CDF), the log of the right tail probability mass) for an
138 * exponential distribution, given value <x>, offset <mu>, and
139 * decay parameter <lambda>.
140 */
141 double
esl_exp_logsurv(double x,double mu,double lambda)142 esl_exp_logsurv(double x, double mu, double lambda)
143 {
144 if (x < mu) return 0.0;
145 return -lambda * (x-mu);
146 }
147
148
149 /* Function: esl_exp_invcdf()
150 *
151 * Purpose: Calculates the inverse of the CDF; given a <cdf> value
152 * $0 <= p < 1$, returns the value $x$ at which the CDF
153 * has that value.
154 */
155 double
esl_exp_invcdf(double p,double mu,double lambda)156 esl_exp_invcdf(double p, double mu, double lambda)
157 {
158 return mu - 1/lambda * log(1. - p);
159 }
160
161
162
163 /* Function: esl_exp_invsurv()
164 *
165 * Purpose: Calculates the inverse of the survivor function, the score
166 * at which the right tail's mass is $0 <= p < 1$, for an
167 * exponential function with parameters <mu> and <lambda>.
168 */
169 double
esl_exp_invsurv(double p,double mu,double lambda)170 esl_exp_invsurv(double p, double mu, double lambda)
171 {
172
173 return mu - 1./lambda * log(p);
174 }
175 /*------------------ end of densities and distributions --------------------*/
176
177
178 /*------------------ end of densities and distributions --------------------*/
179
180
181
182
183 /*****************************************************************
184 * 2. Generic API routines: for general interface w/ histogram module
185 *****************************************************************/
186
187 /* Function: esl_exp_generic_pdf()
188 *
189 * Purpose: Generic-API version of PDF.
190 */
191 double
esl_exp_generic_pdf(double x,void * params)192 esl_exp_generic_pdf(double x, void *params)
193 {
194 double *p = (double *) params;
195 return esl_exp_pdf(x, p[0], p[1]);
196 }
197
198 /* Function: esl_exp_generic_cdf()
199 *
200 * Purpose: Generic-API version of CDF.
201 */
202 double
esl_exp_generic_cdf(double x,void * params)203 esl_exp_generic_cdf(double x, void *params)
204 {
205 double *p = (double *) params;
206 return esl_exp_cdf(x, p[0], p[1]);
207 }
208
209 /* Function: esl_exp_generic_surv()
210 *
211 * Purpose: Generic-API version of survival function.
212 */
213 double
esl_exp_generic_surv(double x,void * params)214 esl_exp_generic_surv(double x, void *params)
215 {
216 double *p = (double *) params;
217 return esl_exp_surv(x, p[0], p[1]);
218 }
219
220 /* Function: esl_exp_generic_invcdf()
221 *
222 * Purpose: Generic-API version of inverse CDF.
223 */
224 double
esl_exp_generic_invcdf(double p,void * params)225 esl_exp_generic_invcdf(double p, void *params)
226 {
227 double *v = (double *) params;
228 return esl_exp_invcdf(p, v[0], v[1]);
229 }
230 /*------------------------- end of generic API --------------------------*/
231
232
233
234 /****************************************************************************
235 * 3. Routines for dumping plots for files
236 ****************************************************************************/
237
238 /* Function: esl_exp_Plot()
239 *
240 * Purpose: Plot some exponential function <func> (for instance,
241 * <esl_exp_pdf()>) for parameters <mu> and <lambda>, for
242 * a range of values x from <xmin> to <xmax> in steps of <xstep>;
243 * output to an open stream <fp> in xmgrace XY input format.
244 *
245 * Returns: <eslOK> on success.
246 *
247 * Throws: <eslEWRITE> on any system write error, such as a filled disk.
248 */
249 int
esl_exp_Plot(FILE * fp,double mu,double lambda,double (* func)(double x,double mu,double lambda),double xmin,double xmax,double xstep)250 esl_exp_Plot(FILE *fp, double mu, double lambda,
251 double (*func)(double x, double mu, double lambda),
252 double xmin, double xmax, double xstep)
253 {
254 double x;
255 for (x = xmin; x <= xmax; x += xstep)
256 if (fprintf(fp, "%f\t%g\n", x, (*func)(x, mu, lambda)) < 0) ESL_EXCEPTION_SYS(eslEWRITE, "exponential plot write failed");
257 if (fprintf(fp, "&\n") < 0) ESL_EXCEPTION_SYS(eslEWRITE, "exponential plot write failed");
258 return eslOK;
259 }
260 /*-------------------- end plot dumping routines ---------------------------*/
261
262
263
264 /****************************************************************************
265 * 4. Routines for sampling
266 ****************************************************************************/
267
268 /* Function: esl_exp_Sample()
269 *
270 * Purpose: Sample an exponential random variate
271 * by the transformation method, given offset <mu>
272 * and decay parameter <lambda>.
273 */
274 double
esl_exp_Sample(ESL_RANDOMNESS * r,double mu,double lambda)275 esl_exp_Sample(ESL_RANDOMNESS *r, double mu, double lambda)
276 {
277 double p, x;
278 p = esl_rnd_UniformPositive(r);
279
280 x = mu - 1./lambda * log(p); // really log(1-p), but if p uniform on 0..1
281 // then so is 1-p.
282 return x;
283 }
284 /*--------------------------- end sampling ---------------------------------*/
285
286
287
288
289 /****************************************************************************
290 * 5. Maximum likelihood fitting
291 ****************************************************************************/
292
293 /* Function: esl_exp_FitComplete()
294 *
295 * Purpose: Given an array of <n> samples <x[0]..x[n-1]>, fit
296 * them to an exponential distribution.
297 * Return maximum likelihood parameters <ret_mu> and <ret_lambda>.
298 *
299 * Args: x - complete exponentially-distributed data [0..n-1]
300 * n - number of samples in <x> (n>0)
301 * ret_mu - lower bound of the distribution (all x_i >= mu)
302 * ret_lambda - RETURN: maximum likelihood estimate of lambda
303 *
304 * Returns: <eslOK> on success.
305 *
306 * Throws: <eslEINVAL> if n=0 (no data).
307 *
308 * Xref: STL9/138.
309 */
310 int
esl_exp_FitComplete(double * x,int n,double * ret_mu,double * ret_lambda)311 esl_exp_FitComplete(double *x, int n, double *ret_mu, double *ret_lambda)
312 {
313 double mu, mean;
314 int i;
315 int status;
316
317 if (!n) ESL_XEXCEPTION(eslEINVAL, "empty data vector provided for exponential fit");
318
319 /* ML mu is the lowest score. mu=x is ok in the exponential. */
320 mu = x[0];
321 for (i = 1; i < n; i++) if (x[i] < mu) mu = x[i];
322
323 mean = 0.;
324 for (i = 0; i < n; i++) mean += x[i] - mu;
325 mean /= (double) n;
326
327 *ret_mu = mu;
328 *ret_lambda = 1./mean; /* ML estimate trivial & analytic */
329 return eslOK;
330
331 ERROR:
332 *ret_mu = 0.0;
333 *ret_lambda = 0.0;
334 return status;
335 }
336
337 /* Function: esl_exp_FitCompleteScale()
338 *
339 * Purpose: Given an array of <n> samples <x[0]..x[n-1]>, fit
340 * them to an exponential distribution of known location
341 * parameter <mu>. Return maximum likelihood scale
342 * parameter <ret_lambda>.
343 *
344 * All $x_i \geq \mu$.
345 *
346 * Args: x - complete exponentially-distributed data [0..n-1]
347 * n - number of samples in <x>
348 * mu - lower bound of the distribution (all x_i >= mu)
349 * ret_lambda - RETURN: maximum likelihood estimate of lambda
350 *
351 * Returns: <eslOK> on success.
352 *
353 * Xref: J1/49.
354 */
355 int
esl_exp_FitCompleteScale(double * x,int n,double mu,double * ret_lambda)356 esl_exp_FitCompleteScale(double *x, int n, double mu, double *ret_lambda)
357 {
358 double mean;
359 int i;
360
361 mean = 0.;
362 for (i = 0; i < n; i++) mean += x[i] - mu;
363 mean /= (double) n;
364
365 *ret_lambda = 1./mean; /* ML estimate trivial & analytic */
366 return eslOK;
367 }
368
369
370 /* Function: esl_exp_FitCompleteBinned()
371 *
372 * Purpose: Fit a complete exponential distribution to the observed
373 * binned data in a histogram <g>, where each
374 * bin i holds some number of observed samples x with values from
375 * lower bound l to upper bound u (that is, $l < x \leq u$);
376 * find maximum likelihood parameters $\mu,\lambda$ and
377 * return them in <*ret_mu>, <*ret_lambda>.
378 *
379 * If the binned data in <g> were set to focus on
380 * a tail by virtual censoring, the "complete" exponential is
381 * fitted to this tail. The caller then also needs to
382 * remember what fraction of the probability mass was in this
383 * tail.
384 *
385 * The ML estimate for $mu$ is the smallest observed
386 * sample. For complete data, <ret_mu> is generally set to
387 * the smallest observed sample value, except in the
388 * special case of a "rounded" complete dataset, where
389 * <ret_mu> is set to the lower bound of the smallest
390 * occupied bin. For tails, <ret_mu> is set to the cutoff
391 * threshold <phi>, where we are guaranteed that <phi> is
392 * at the lower bound of a bin (by how the histogram
393 * object sets tails).
394 *
395 * The ML estimate for <ret_lambda> has an analytical
396 * solution, so this routine is fast.
397 *
398 * If all the data are in one bin, the ML estimate of
399 * $\lambda$ will be $\infty$. This is mathematically correct,
400 * but is probably a situation the caller wants to avoid, perhaps
401 * by choosing smaller bins.
402 *
403 * This function currently cannot fit an exponential tail
404 * to truly censored, binned data, because it assumes that
405 * all bins have equal width, but in true censored data, the
406 * lower cutoff <phi> may fall anywhere in the first bin.
407 *
408 * Returns: <eslOK> on success.
409 *
410 * Throws: <eslEINVAL> if dataset is true-censored.
411 */
412 int
esl_exp_FitCompleteBinned(ESL_HISTOGRAM * g,double * ret_mu,double * ret_lambda)413 esl_exp_FitCompleteBinned(ESL_HISTOGRAM *g, double *ret_mu, double *ret_lambda)
414 {
415 int i;
416 double ai, bi, delta;
417 double sa, sb;
418 double mu = 0.;
419
420 if (g->dataset_is == COMPLETE)
421 {
422 if (g->is_rounded) mu = esl_histogram_Bin2LBound(g, g->imin);
423 else mu = g->xmin;
424 }
425 else if (g->dataset_is == VIRTUAL_CENSORED) /* i.e., we'll fit to tail */
426 mu = g->phi;
427 else if (g->dataset_is == TRUE_CENSORED)
428 ESL_EXCEPTION(eslEINVAL, "can't fit true censored dataset");
429
430 delta = g->w;
431 sa = sb = 0.;
432 for (i = g->cmin; i <= g->imax; i++) /* for each occupied bin */
433 {
434 if (g->obs[i] == 0) continue;
435 ai = esl_histogram_Bin2LBound(g,i);
436 bi = esl_histogram_Bin2UBound(g,i);
437 sa += g->obs[i] * (ai-mu);
438 sb += g->obs[i] * (bi-mu);
439 }
440 *ret_mu = mu;
441 *ret_lambda = 1/delta * (log(sb) - log(sa));
442 return eslOK;
443 }
444
445
446 /****************************************************************************
447 * 6. Stats driver
448 ****************************************************************************/
449 #ifdef eslEXPONENTIAL_STATS
450 /* Compiles statistics on the accuracy of ML estimation of an exponential tail.
451 * compile: gcc -g -O2 -Wall -I. -L. -o stats -DeslEXPONENTIAL_STATS esl_exponential.c -leasel -lm
452 * run: ./stats > stats.out
453 *
454 * Output is, for each trial:
455 * <trial #> <fitted mu> <fitted lambda>
456 *
457 * To get mean, stddev of lambda estimates:
458 * % ./stats | avg -f2
459 */
460 #include <stdio.h>
461
462 #include "easel.h"
463 #include "esl_getopts.h"
464 #include "esl_random.h"
465 #include "esl_exponential.h"
466
467 int
main(int argc,char ** argv)468 main(int argc, char **argv)
469 {
470 ESL_RANDOMNESS *r = esl_randomness_Create(0);
471 int ntrials; /* number of estimates to gather */
472 int N; /* number of samples collected to make each estimate */
473 double mu, lambda; /* parametric location, scale */
474 double est_mu, est_lambda; /* estimated location, scale */
475 int trial;
476 int i;
477 double *x;
478
479 /* Configuration: (change & recompile as needed)
480 */
481 ntrials = 1000;
482 mu = 0.;
483 lambda = 0.693;
484 N = 95;
485
486 x = malloc(sizeof(double) *N);
487 for (trial = 0; trial < ntrials; trial++)
488 {
489 for (i = 0; i < N; i++)
490 x[i] = esl_exp_Sample(r, mu, lambda);
491 esl_exp_FitComplete(x, N, &est_mu, &est_lambda);
492
493 /*
494 est_mu = mu;
495 esl_exp_FitCompleteScale(x, N, est_mu, &est_lambda);
496 */
497 printf("%4d %8.4f %8.4f\n", i, est_mu, est_lambda);
498 }
499 free(x);
500 return 0;
501 }
502 #endif /*eslEXPONENTIAL_STATS*/
503
504
505
506
507
508 /****************************************************************************
509 * 7. Unit tests
510 ****************************************************************************/
511
512 /****************************************************************************
513 * 8. Test driver
514 ****************************************************************************/
515 #ifdef eslEXPONENTIAL_TESTDRIVE
516 /* Compile:
517 gcc -g -Wall -I. -L. -o test -DeslEXPONENTIAL_TESTDRIVE esl_exponential.c -leasel -lm
518 */
519 #include <stdio.h>
520 #include <stdlib.h>
521 #include <string.h>
522
523 #include "easel.h"
524 #include "esl_random.h"
525 #include "esl_histogram.h"
526 #include "esl_exponential.h"
527
528 int
main(int argc,char ** argv)529 main(int argc, char **argv)
530 {
531 ESL_HISTOGRAM *h;
532 ESL_RANDOMNESS *r;
533 double mu = 10.0;
534 double lambda = 1.0;
535 int n = 10000;
536 double binwidth = 0.1;
537 double emu, elambda;
538 int i;
539 double x;
540 double *data;
541 int ndata;
542
543 int opti;
544 int be_verbose = FALSE;
545 char *plotfile = NULL;
546 FILE *pfp = stdout;
547 int plot_pdf = FALSE;
548 int plot_logpdf = FALSE;
549 int plot_cdf = FALSE;
550 int plot_logcdf = FALSE;
551 int plot_surv = FALSE;
552 int plot_logsurv = FALSE;
553 int xmin_set = FALSE;
554 double xmin;
555 int xmax_set = FALSE;
556 double xmax;
557 int xstep_set = FALSE;
558 double xstep;
559
560 for (opti = 1; opti < argc && *(argv[opti]) == '-'; opti++)
561 {
562 if (strcmp(argv[opti], "-m") == 0) mu = atof(argv[++opti]);
563 else if (strcmp(argv[opti], "-l") == 0) lambda = atof(argv[++opti]);
564 else if (strcmp(argv[opti], "-n") == 0) n = atoi(argv[++opti]);
565 else if (strcmp(argv[opti], "-o") == 0) plotfile = argv[++opti];
566 else if (strcmp(argv[opti], "-v") == 0) be_verbose = TRUE;
567 else if (strcmp(argv[opti], "-w") == 0) binwidth = atof(argv[++opti]);
568 else if (strcmp(argv[opti], "-C") == 0) plot_cdf = TRUE;
569 else if (strcmp(argv[opti], "-LC") == 0) plot_logcdf = TRUE;
570 else if (strcmp(argv[opti], "-P") == 0) plot_pdf = TRUE;
571 else if (strcmp(argv[opti], "-LP") == 0) plot_logpdf = TRUE;
572 else if (strcmp(argv[opti], "-S") == 0) plot_surv = TRUE;
573 else if (strcmp(argv[opti], "-LS") == 0) plot_logsurv = TRUE;
574 else if (strcmp(argv[opti], "-XL") == 0) { xmin_set = TRUE; xmin = atof(argv[++opti]); }
575 else if (strcmp(argv[opti], "-XH") == 0) { xmax_set = TRUE; xmax = atof(argv[++opti]); }
576 else if (strcmp(argv[opti], "-XS") == 0) { xstep_set = TRUE; xstep = atof(argv[++opti]); }
577 else esl_fatal("bad option");
578 }
579
580 if (be_verbose)
581 printf("Parametric: mu = %f lambda = %f\n", mu, lambda);
582
583 r = esl_randomness_Create(0);
584 h = esl_histogram_CreateFull(mu, 100., binwidth);
585 if (plotfile != NULL) {
586 if ((pfp = fopen(plotfile, "w")) == NULL) esl_fatal("Failed to open plotfile");
587 }
588 if (! xmin_set) xmin = mu;
589 if (! xmax_set) xmax = mu+20* (1./lambda);
590 if (! xstep_set) xstep = 0.1;
591
592 for (i = 0; i < n; i++)
593 {
594 x = esl_exp_Sample(r, mu, lambda);
595 esl_histogram_Add(h, x);
596 }
597 esl_histogram_GetData(h, &data, &ndata);
598
599 esl_exp_FitComplete(data, ndata, &emu, &elambda);
600 if (be_verbose) printf("Complete data fit: mu = %f lambda = %f\n", emu, elambda);
601 if (fabs( (emu-mu)/mu ) > 0.01) esl_fatal("Error in (complete) fitted mu > 1%\n");
602 if (fabs( (elambda-lambda)/lambda ) > 0.10) esl_fatal("Error in (complete) fitted lambda > 10%\n");
603
604 esl_exp_FitCompleteBinned(h, &emu, &elambda);
605 if (be_verbose) printf("Binned data fit: mu = %f lambda = %f\n", emu, elambda);
606 if (fabs( (emu-mu)/mu ) > 0.01) esl_fatal("Error in (binned) fitted mu > 1%\n");
607 if (fabs( (elambda-lambda)/lambda ) > 0.10) esl_fatal("Error in (binned) fitted lambda > 10%\n");
608
609 if (plot_pdf) esl_exp_Plot(pfp, mu, lambda, &esl_exp_pdf, xmin, xmax, xstep);
610 if (plot_logpdf) esl_exp_Plot(pfp, mu, lambda, &esl_exp_logpdf, xmin, xmax, xstep);
611 if (plot_cdf) esl_exp_Plot(pfp, mu, lambda, &esl_exp_cdf, xmin, xmax, xstep);
612 if (plot_logcdf) esl_exp_Plot(pfp, mu, lambda, &esl_exp_logcdf, xmin, xmax, xstep);
613 if (plot_surv) esl_exp_Plot(pfp, mu, lambda, &esl_exp_surv, xmin, xmax, xstep);
614 if (plot_logsurv) esl_exp_Plot(pfp, mu, lambda, &esl_exp_logsurv, xmin, xmax, xstep);
615
616 if (plotfile != NULL) fclose(pfp);
617
618 esl_randomness_Destroy(r);
619 esl_histogram_Destroy(h);
620 return 0;
621 }
622 #endif /*eslEXPONENTIAL_TESTDRIVE*/
623
624
625 /****************************************************************************
626 * 9. Example
627 ****************************************************************************/
628 #ifdef eslEXPONENTIAL_EXAMPLE
629 /*::cexcerpt::exp_example::begin::*/
630 #include <stdio.h>
631 #include "easel.h"
632 #include "esl_random.h"
633 #include "esl_histogram.h"
634 #include "esl_exponential.h"
635
636 int
main(int argc,char ** argv)637 main(int argc, char **argv)
638 {
639 double mu = -50.0;
640 double lambda = 0.5;
641 ESL_RANDOMNESS *r = esl_randomness_Create(0);
642 ESL_HISTOGRAM *h = esl_histogram_CreateFull(mu, 100., 0.1);
643 int n = 10000;
644 double emu, elambda;
645 int i;
646 double x;
647 double *data;
648 int ndata;
649
650 for (i = 0; i < n; i++)
651 {
652 x = esl_exp_Sample(r, mu, lambda);
653 esl_histogram_Add(h, x);
654 }
655 esl_histogram_GetData(h, &data, &ndata);
656
657 /* Plot the empirical (sampled) and expected survivals */
658 esl_histogram_PlotSurvival(stdout, h);
659 esl_exp_Plot(stdout, mu, lambda,
660 &esl_exp_surv, h->xmin, h->xmax, 0.1);
661
662 /* ML fit to complete data, and plot fitted survival curve */
663 esl_exp_FitComplete(data, ndata, &emu, &elambda);
664 esl_exp_Plot(stdout, emu, elambda,
665 &esl_exp_surv, h->xmin, h->xmax, 0.1);
666
667 /* ML fit to binned data, plot fitted survival curve */
668 esl_exp_FitCompleteBinned(h, &emu, &elambda);
669 esl_exp_Plot(stdout, emu, elambda,
670 &esl_exp_surv, h->xmin, h->xmax, 0.1);
671
672 esl_randomness_Destroy(r);
673 esl_histogram_Destroy(h);
674 return 0;
675 }
676 /*::cexcerpt::exp_example::end::*/
677 #endif /*eslEXPONENTIAL_EXAMPLE*/
678
679