1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: test_StdDistr.c *
8 * *
9 * Compare CDF, PDF and derivatives of PDF of varios distributions *
10 * with Mathematica(TM) output. *
11 * *
12 *****************************************************************************
13 * *
14 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
15 * Department of Statistics and Mathematics, WU Wien, Austria *
16 * *
17 * This program is free software; you can redistribute it and/or modify *
18 * it under the terms of the GNU General Public License as published by *
19 * the Free Software Foundation; either version 2 of the License, or *
20 * (at your option) any later version. *
21 * *
22 * This program is distributed in the hope that it will be useful, *
23 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
24 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
25 * GNU General Public License for more details. *
26 * *
27 * You should have received a copy of the GNU General Public License *
28 * along with this program; if not, write to the *
29 * Free Software Foundation, Inc., *
30 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
31 * *
32 *****************************************************************************/
33
34 /*---------------------------------------------------------------------------*/
35
36 /* #define DEBUG 1 */
37
38 /*---------------------------------------------------------------------------*/
39 #include <ctype.h>
40 #include <stdio.h>
41 #include <stdlib.h>
42 #include <string.h>
43 #include <time.h>
44
45 #include <config.h>
46 #include <unuran.h>
47
48 /*---------------------------------------------------------------------------*/
49
50 /* PDF test: maximal difference allowed (relative to maximal value) */
51 #define MAX_REL_DIFF 1.e-10
52
53 /* mode test: */
54 #define EPSX 1.e-5 /* stepsize */
55 #define EPSY 1.e-12 /* maximal allowed relative error */
56
57 /*---------------------------------------------------------------------------*/
58
59 /* name of data file */
60 static const char datafile[] = "t_StdDistr.data";
61
62 /* file handles */
63 static FILE *DATA; /* file with data */
64 static FILE *UNURANLOG; /* unuran log file */
65 static FILE *TESTLOG; /* test log file */
66
67 /*---------------------------------------------------------------------------*/
68
69 static int _unur_test_StdDistr (int *n_failed);
70 /* main loop for testing standard distributions */
71
72 static int test_cdf_pdf( UNUR_DISTR *distr, char *datafile );
73 /* compare CDF, PDF and dPDF for continuous univariate distributions */
74
75 static int modetest_cont( UNUR_DISTR *distr);
76 static int modetest_discr( UNUR_DISTR *distr);
77 /* test mode of distribution */
78
79 /*---------------------------------------------------------------------------*/
80
81 #define MAX(a,b) ( ((a) < (b)) ? (b) : (a))
82 #define MIN(a,b) ( ((a) > (b)) ? (b) : (a))
83
84 /*---------------------------------------------------------------------------*/
85
main()86 int main()
87 {
88 /* number of failed tests */
89 int n_failed = 0;
90
91 /* open log file for unuran and set output stream for unuran messages */
92 if ( (UNURANLOG = fopen( "t_StdDistr_unuran.log","w" )) == NULL )
93 exit (EXIT_FAILURE);
94 unur_set_stream( UNURANLOG );
95
96 /* open log file for testing */
97 if ( (TESTLOG = fopen( "t_StdDistr_test.log","w" )) == NULL )
98 exit (EXIT_FAILURE);
99
100 /* write header into log file */
101 {
102 time_t started;
103 fprintf(TESTLOG,"\nUNU.RAN - Universal Non-Uniform RANdom number generator\n\n");
104 if (time( &started ) != -1)
105 fprintf(TESTLOG,"%s",ctime(&started));
106 fprintf(TESTLOG,"\n=======================================================\n\n");
107 fprintf(TESTLOG,"(Search for string \"distr=\" to find new section.)\n\n");
108 }
109
110 /* open data file */
111 if ( (DATA = fopen( datafile,"r" )) == NULL ) {
112 printf("ERROR: could not open file %s \n", datafile);
113 fprintf(TESTLOG,"ERROR: could not open file %s \n", datafile);
114 exit (77); /* ignore this error */
115 }
116
117 /* run tests on all distributions */
118 while (_unur_test_StdDistr(&n_failed)==TRUE);
119
120 /* close files */
121 fclose (UNURANLOG);
122 fclose (TESTLOG);
123 fclose (DATA);
124
125 /* end */
126 if (n_failed > 0) {
127 printf("[StdDistr --> failed]\n");
128 exit (EXIT_FAILURE);
129 }
130 else {
131 printf("[StdDistr --> ok]\n");
132 exit (EXIT_SUCCESS);
133 }
134 } /* end of main() */
135
136 /*---------------------------------------------------------------------------*/
137
138 int
_unur_test_StdDistr(int * n_failed)139 _unur_test_StdDistr (int *n_failed)
140 {
141 #define BUFSIZE 1024 /* size of line buffer */
142
143 char buffer[BUFSIZE]; /* buffer for reading line */
144 char dstr[BUFSIZE]; /* distribution string */
145
146 UNUR_DISTR *distr; /* distribution object */
147
148 /* find next function string */
149 while (1) {
150 /* read next line */
151 fgets(buffer, BUFSIZE, DATA);
152 if (feof(DATA)) return FALSE;
153 if (strncmp(buffer,"distr=",6)==0) break;
154 }
155
156 /* get function string */
157 strcpy( dstr, buffer+6 );
158 /* remove newline character */
159 dstr[strlen(dstr)-1] = '\0';
160
161 /* make distribution object with given function as PDF */
162 if ( (distr = unur_str2distr(dstr)) == NULL ) {
163 printf("ERROR: syntax error in \"%s\"\n", dstr);
164 fprintf(TESTLOG,"ERROR: syntax error in \"%s\"\n", dstr);
165 exit (EXIT_FAILURE);
166 }
167
168 /* now run test */
169 if (test_cdf_pdf( distr,dstr ) == UNUR_FAILURE)
170 (*n_failed)++;
171
172 /* free memory */
173 unur_distr_free(distr);
174
175 return TRUE;
176
177 #undef BUFSIZE
178 } /* end of _unur_test_StdDistr() */
179
180 /*---------------------------------------------------------------------------*/
181
182 int
test_cdf_pdf(UNUR_DISTR * distr,char * distrAPI)183 test_cdf_pdf( UNUR_DISTR *distr, char *distrAPI )
184 /*----------------------------------------------------------------------*/
185 /* test CDF, PDF and derivative of PDF by comparing to data */
186 /* read from file created by Mathematica. */
187 /* */
188 /* the ordering in this file is: */
189 /* n par[1] par[2] ... par[n] x CDF PDF dPDF */
190 /* */
191 /* parameters: */
192 /* distr ... pointer to distribution object */
193 /* distrAPI ... string that contains name for distribution */
194 /* */
195 /* return: */
196 /* UNUR_SUCCESS ... if test was successful */
197 /* UNUR_FAILURE ... failed (difference too large) */
198 /* UNUR_ERR_NULL ... could not run test */
199 /*----------------------------------------------------------------------*/
200 {
201 #define BUFSIZE 1024 /* size of line buffer */
202 #define MAX_FPARAMS 10 /* maximal number of parameters for distribution */
203
204 char buffer[BUFSIZE]; /* buffer for reading line */
205 char *ptr_buffer; /* pointer into buffer */
206
207 int n_fparams; /* number of parameters for distribution */
208 double fparams[10]; /* array for parameters for distribution */
209 double x; /* double argument */
210 int k = 0; /* integer argument */
211 double CDF_e, PDF_e, dPDF_e; /* expected values for CDF, PDF and derivative of PDF at x */
212 double CDF_o, PDF_o, dPDF_o; /* observed values for CDF, PDF and derivative of PDF at x */
213 double CDF_d, PDF_d, dPDF_d; /* differences between observed and expected values */
214
215 double CDF_md, PDF_md, dPDF_md; /* maximal difference (absolute) */
216 double CDF_me, PDF_me, dPDF_me; /* maximal expected absolute values */
217
218 int have_CDF, have_PDF, have_dPDF, have_upd_pdfarea;
219
220 int is_DISCR; /* 1 if discrete distribution, 0 otherwise */
221 int n_failed = 0; /* number of failed tests */
222
223 const char *dname; /* name of distribution */
224
225 int i;
226
227 /* discrete distribution ? */
228 is_DISCR = unur_distr_is_discr( distr );
229
230 /* name of distribution */
231 dname = unur_distr_get_name(distr);
232
233 /* initialize */
234 CDF_md = PDF_md = dPDF_md = 0.;
235 CDF_me = PDF_me = dPDF_me = 0.;
236
237 /* check existence of CDF, PDF and dPDF */
238 if (is_DISCR==TRUE) {
239 have_CDF = (unur_distr_discr_get_cdf(distr)) ? TRUE : FALSE;
240 have_PDF = (unur_distr_discr_get_pmf(distr)) ? TRUE : FALSE;
241 have_dPDF = FALSE;
242 }
243 else { /* is_CONT */
244 have_CDF = (unur_distr_cont_get_cdf(distr)) ? TRUE : FALSE;
245 have_PDF = (unur_distr_cont_get_pdf(distr)) ? TRUE : FALSE;
246 have_dPDF = (unur_distr_cont_get_dpdf(distr)) ? TRUE : FALSE;
247 }
248
249 /* check whether unur_distr_cont_upd_pdfarea() works */
250 if (is_DISCR==TRUE)
251 have_upd_pdfarea = (unur_distr_discr_upd_pmfsum(distr)==UNUR_SUCCESS) ? TRUE : FALSE;
252 else /* is_CONT */
253 have_upd_pdfarea = (unur_distr_cont_upd_pdfarea(distr)==UNUR_SUCCESS) ? TRUE : FALSE;
254
255 if (have_upd_pdfarea==FALSE) {
256 /* if we cannot update the area below the PDF, then the
257 given PDF is not a "real" PDF, i.e. not normalized */
258 have_PDF = FALSE;
259 have_dPDF = FALSE;
260 }
261
262 /* print info into log file */
263 fprintf(TESTLOG,"%s: distr= %s ...\n", dname, distrAPI);
264
265 if (have_CDF==FALSE)
266 fprintf(TESTLOG,"%s: no CDF!\n", dname);
267 if (have_PDF==FALSE)
268 fprintf(TESTLOG,"%s: no PDF!\n", dname);
269 if (have_dPDF==FALSE)
270 fprintf(TESTLOG,"%s: no dPDF!\n", dname);
271
272 /* read data file */
273 while (1) {
274
275 /* read next line */
276 fgets(buffer, BUFSIZE-1, DATA);
277 if (feof(DATA) || isspace(buffer[0])) break;
278
279 ptr_buffer = buffer;
280
281 /* get number of parameters */
282 n_fparams = strtol(buffer, &ptr_buffer, 10);
283 if (n_fparams < 0 || n_fparams >= MAX_FPARAMS) {
284 printf("%s: ERROR: invalid number of parameters for distribution: %d \n", dname,n_fparams);
285 fprintf(TESTLOG,"%s: ERROR: invalid number of parameters for distribution: %d \n", dname,n_fparams);
286 return UNUR_ERR_NULL;
287 }
288
289 /* read parameters */
290 for (i=0; i<n_fparams; i++) {
291 fparams[i] = strtod( ptr_buffer, &ptr_buffer );
292 }
293
294 /* read argument */
295 x = strtod( ptr_buffer, &ptr_buffer );
296 if (is_DISCR)
297 /* integer required, round */
298 k = (int)(x+0.5);
299
300 /* read CDF */
301 CDF_e = strtod( ptr_buffer, &ptr_buffer );
302
303 /* read PDF */
304 PDF_e = strtod( ptr_buffer, &ptr_buffer );
305
306 /* read dPDF */
307 dPDF_e = strtod( ptr_buffer, &ptr_buffer );
308
309 #ifdef DEBUG
310 /* print data */
311 fprintf(TESTLOG, "%s: n_fparams = %d\n", dname, n_fparams);
312
313 /* print parameters */
314 for (i=0; i<n_fparams; i++)
315 fprintf(TESTLOG, "%s: \t%d: %g\n", dname, i, fparams[i]);
316
317 /* print argument x (or k) */
318 if (is_DISCR)
319 fprintf(TESTLOG, "%s: expected k = %d:\t", dname, k);
320 else /* is_CONT */
321 fprintf(TESTLOG, "%s: expected x = %g:\t", dname, x);
322
323 /* print CDF, PDF and derivative of PDF at x */
324 fprintf(TESTLOG, "%g, %g, %g\n", CDF_e, PDF_e, dPDF_e);
325 #endif
326
327 if (is_DISCR) {
328 /* set parameters for distribution */
329 unur_distr_discr_set_pmfparams(distr,fparams,n_fparams);
330 if (have_upd_pdfarea) unur_distr_discr_upd_pmfsum(distr);
331
332 /* compute CDF, PDF and derivative of PDF at x */
333 CDF_o = (have_CDF) ? unur_distr_discr_eval_cdf(x, distr) : 0.;
334 PDF_o = (have_PDF) ? unur_distr_discr_eval_pmf(x, distr) : 0.;
335 dPDF_o = 0.;
336 }
337
338 else { /* is_CONT */
339 /* set parameters for distribution */
340 unur_distr_cont_set_pdfparams(distr,fparams,n_fparams);
341 if (have_upd_pdfarea) unur_distr_cont_upd_pdfarea(distr);
342
343 /* compute CDF, PDF and derivative of PDF at x */
344 CDF_o = (have_CDF) ? unur_distr_cont_eval_cdf(x, distr) : 0.;
345 PDF_o = (have_PDF) ? unur_distr_cont_eval_pdf(x, distr) : 0.;
346 dPDF_o = (have_dPDF) ? unur_distr_cont_eval_dpdf(x, distr) : 0.;
347 }
348
349 /* compute differnces */
350 CDF_d = CDF_o - CDF_e;
351 PDF_d = PDF_o - PDF_e;
352 dPDF_d = dPDF_o - dPDF_e;
353
354 #ifdef DEBUG
355 /* print argument x (or k) */
356 if (is_DISCR)
357 fprintf(TESTLOG, "%s: observed k = %d:\t", dname, k);
358 else /* is_CONT */
359 fprintf(TESTLOG, "%s: observed x = %g:\t", dname, x);
360
361 /* print CDF, PDF and derivative of PDF at x */
362 fprintf(TESTLOG, "%g, %g, %g\n",CDF_o,PDF_o,dPDF_o);
363
364 /* print differences */
365 fprintf(TESTLOG, "%s: diff x = %g:\t", dname,x);
366 fprintf(TESTLOG, "%g, %g, %g\n",CDF_d, PDF_d, dPDF_d);
367 fprintf(TESTLOG, "%s:\n", dname);
368 #endif
369
370 /* absolute values */
371 CDF_d = fabs(CDF_d);
372 PDF_d = fabs(PDF_d);
373 dPDF_d = fabs(dPDF_d);
374
375 CDF_e = fabs(CDF_e);
376 PDF_e = fabs(PDF_e);
377 dPDF_e = fabs(dPDF_e);
378
379 /* maximal differences */
380 if (CDF_d > CDF_md) CDF_md = CDF_d;
381 if (PDF_d > PDF_md) PDF_md = PDF_d;
382 if (dPDF_d > dPDF_md) dPDF_md = dPDF_d;
383
384 /* maximal values */
385 if (CDF_e > CDF_me) CDF_me = CDF_e;
386 if (PDF_e > PDF_me) PDF_me = PDF_e;
387 if (dPDF_e > dPDF_me) dPDF_me = dPDF_e;
388
389 /* test mode of distribution */
390 if (is_DISCR) {
391 if (modetest_discr(distr) == UNUR_FAILURE)
392 ++n_failed;
393 }
394 else { /* is_CONT */
395 if (modetest_cont(distr) == UNUR_FAILURE)
396 ++n_failed;
397 }
398
399 }
400
401 /* print info on screen */
402 fprintf(TESTLOG, "%s: \tmaximal difference:\n", dname);
403
404 if (have_CDF) {
405 fprintf(TESTLOG,"%s: \t\tCDF = %g ... ", dname,CDF_md);
406 if (CDF_md > MAX_REL_DIFF * CDF_me) {
407 fprintf(TESTLOG, "failed!!\n");
408 ++n_failed;
409 }
410 else
411 fprintf(TESTLOG, "ok\n");
412 }
413
414 if (have_PDF) {
415 fprintf(TESTLOG,"%s: \t\tPDF = %g (rel = %g) ... ", dname,PDF_md,PDF_md/PDF_me);
416 if (PDF_md > MAX_REL_DIFF * PDF_me) {
417 fprintf(TESTLOG, "failed!!\n");
418 ++n_failed;
419 }
420 else
421 fprintf(TESTLOG, "ok\n");
422 }
423
424 if (have_dPDF) {
425 fprintf(TESTLOG,"%s: \t\tdPDF = %g (rel = %g) ... ", dname,dPDF_md,dPDF_md/dPDF_me);
426 if (dPDF_md > MAX_REL_DIFF * dPDF_me) {
427 fprintf(TESTLOG, "failed!!\n");
428 ++n_failed;
429 }
430 else
431 fprintf(TESTLOG, "ok\n");
432 }
433
434 fprintf(TESTLOG, "\n----------------------------------------\n\n");
435
436 /* print result on screen */
437 printf("%-17s ... ", dname );
438 if (n_failed > 0) {
439 printf("failed!!\n");
440 return UNUR_FAILURE;
441 }
442 else {
443 printf("ok\n");
444 return UNUR_SUCCESS;
445 }
446
447 #undef BUFSIZE
448 #undef MAX_FPARAMS
449 } /* end of test_cdf_pdf() */
450
451 /*---------------------------------------------------------------------------*/
452
modetest_cont(UNUR_DISTR * distr)453 int modetest_cont( UNUR_DISTR *distr)
454 /*----------------------------------------------------------------------*/
455 /* Tests whether the mode of continuous distribution is correct. */
456 /* */
457 /* parameters: */
458 /* distr ... pointer to distribution object */
459 /* */
460 /* return: */
461 /* UNUR_SUCCESS ... if test was successful */
462 /* UNUR_FAILURE ... failed (difference too large) */
463 /* UNUR_ERR_NULL ... could not run test */
464 /*----------------------------------------------------------------------*/
465 {
466 double m, fm; /* mode and value of PDF at mode */
467 double x, fx; /* argument x and value of PDF at x */
468 double domain[2]; /* boundaries of the domain of distribution */
469 const char *dname; /* name of distribution */
470 int n_failed = 0; /* number of failed tests */
471 int i; /* loop variable */
472
473 /* name of distribution */
474 dname = unur_distr_get_name(distr);
475
476 /* get mode */
477 if (unur_distr_cont_upd_mode(distr)!=UNUR_SUCCESS ) {
478 /* cannot read mode */
479 printf("%s: ERROR: cannot get mode\n", dname);
480 fprintf(TESTLOG,"%s: ERROR: cannot get mode\n", dname);
481 return UNUR_ERR_NULL;
482 }
483 m = unur_distr_cont_get_mode(distr);
484
485 /* Correct possible problems if the mode m is on the boundary of domain */
486 unur_distr_cont_get_domain(distr,domain, domain+1);
487 if(domain[1] - m < 1.e-11)
488 m = MIN( m-EPSX/2., m*(m<0?(1.+EPSX/2.):1.-EPSX/2.) );
489 if(m - domain[0] < 1.e-11)
490 m = MAX( m+EPSX/2., m*(m>0?(1.+EPSX/2.):1.-EPSX/2.) );
491
492 /* evaluate PDF at mode */
493 fm = unur_distr_cont_eval_pdf(m, distr);
494
495 #ifdef DEBUG
496 fprintf(TESTLOG,"%s: mode: m = %.20g, f(m) = %.20g\n", dname,m,fm);
497 #endif
498
499 /* test PDF left and right of mode */
500 for (i=0; i<2; i++) {
501
502 if (i==0)
503 /* right of mode */
504 x = MAX( m+EPSX, m*(m>0?(1.+EPSX):1.-EPSX) );
505
506 else
507 /* left of mode */
508 x = MIN( m-EPSX, m*(m<0?(1.+EPSX):1.-EPSX) );
509
510 /* evaluate PDF */
511 fx = unur_distr_cont_eval_pdf(x, distr);
512
513 #ifdef DEBUG
514 fprintf(TESTLOG,"%s: x = %.20g, f(x) = %.20g", dname,x,fx);
515 #endif
516
517 if(fm * (1.+EPSY) < fx) {
518 #ifdef DEBUG
519 fprintf(TESTLOG," ... failed! f(mode) not maximal!");
520 #endif
521 ++n_failed;
522 }
523
524 #ifdef DEBUG
525 fprintf(TESTLOG,"\n");
526 #endif
527 }
528
529 #ifdef DEBUG
530 fprintf(TESTLOG,"%s:\n",dname);
531 #endif
532
533 /* end */
534 return ((n_failed > 0) ? UNUR_FAILURE : UNUR_SUCCESS);
535
536 } /* end of modetest_cont() */
537
538 /*---------------------------------------------------------------------------*/
539
modetest_discr(UNUR_DISTR * distr)540 int modetest_discr( UNUR_DISTR *distr)
541 /*----------------------------------------------------------------------*/
542 /* Tests whether the mode of discr distribution is correct. */
543 /* */
544 /* parameters: */
545 /* distr ... pointer to distribution object */
546 /* */
547 /* return: */
548 /* UNUR_SUCCESS ... if test was successful */
549 /* UNUR_FAILURE ... failed (difference too large) */
550 /* UNUR_ERR_NULL ... could not run test */
551 /*----------------------------------------------------------------------*/
552 {
553 int m, x; /* mode and other argument for PDF */
554 double fm, fx; /* value of PDF at m and x */
555 const char *dname; /* name of distribution */
556 int n_failed = 0; /* number of failed tests */
557 int i; /* loop variable */
558
559 /* name of distribution */
560 dname = unur_distr_get_name(distr);
561
562 /* get mode */
563 if (unur_distr_discr_upd_mode(distr)!=UNUR_SUCCESS ) {
564 /* cannot read mode */
565 printf("%s: ERROR: cannot get mode\n", dname);
566 fprintf(TESTLOG,"%s: ERROR: cannot get mode\n", dname);
567 return UNUR_ERR_NULL;
568 }
569 m = unur_distr_discr_get_mode(distr);
570
571 /* evaluate PMF at mode */
572 fm = unur_distr_discr_eval_pmf(m, distr);
573
574 #ifdef DEBUG
575 fprintf(TESTLOG,"%s: mode: m = %d, f(m) = %.20g\n", dname,m,fm);
576 #endif
577
578 /* test PMF left and right of mode */
579 for (i=-1; i<2; i+=2 ) {
580
581 /* left and right of mode */
582 x = m + i;
583
584 /* evaluate PMF */
585 fx = unur_distr_discr_eval_pmf(x, distr);
586
587 #ifdef DEBUG
588 fprintf(TESTLOG,"%s: x = %d, f(x) = %.20g", dname,x,fx);
589 #endif
590
591 if(fm * (1.+EPSY) < fx) {
592 #ifdef DEBUG
593 fprintf(TESTLOG," ... failed! f(mode) not maximal!");
594 #endif
595 ++n_failed;
596 }
597
598 #ifdef DEBUG
599 fprintf(TESTLOG,"\n");
600 #endif
601 }
602
603 #ifdef DEBUG
604 fprintf(TESTLOG,"%s:\n",dname);
605 #endif
606
607 /* end */
608 return ((n_failed > 0) ? UNUR_FAILURE : UNUR_SUCCESS);
609
610 } /* end of modetest_cont() */
611
612 /*---------------------------------------------------------------------------*/
613
614