1 1 SMOOTH: COMPUTING MONOTONE REGRESSION SPLINES FOR SMOOTHED BOOTSTRAPING CLARK GAYLORD AND DONALD E. RAMIREZ DEPARTMENT OF MATHEMATICS MATH-ASTRO BUILDING UNIVERSITY OF VIRGINIA CHARLOTTESVILLE, VIRGINIA 22903 TECHNICAL REPORT (JUNE 29, 1989) 1 INTRODUCTION THIS TECHNICAL REPORT PRESENTS THE PROGRAM SMOOTH WHICH IS DEVELOPED IN OUR PAPER: MONOTONE REGRESSION SPLINES FOR SMOOTHED BOOTSTRAPING SMOOTH IS A FORTRAN 77 PROGRAM WHICH COMPUTES A REGRESSION SPLINE USING A WEIGHTED INNER PRODUCT TO ACHIEVE MONOTONICITY. THE PROGRAM SMOOTH REQUIRES THAT THE IMSL LIBRARY BE AVAILABLE. INDEX APPENDIX A CONTAINS THE RAW DATA SET WITH N=100 APPENDIX B CONTAINS THE TERMINAL OUTPUT FROM SMOOTH APPENDIX C CONTAINS THE FULL OUTPUT FROM THE PROGRAM SMOOTH APPENDIX D CONTAINS THE FORTRAN PROGRAM SMOOTH APPENDIX E CONTAINS THE CALL STATEMENTS FROM SMOOTH APPENDIX A: THE ORIGINAL DATA 0.57837 0.81219 0.44829 0.37628 0.25938 0.58797 0.44336 0.61926 0.69772 0.28014 0.66886 0.61229 0.29274 0.23483 0.36823 0.30109 0.58663 0.65648 0.11257 0.97886 0.38133 0.60836 0.64048 0.40594 0.31956 0.42196 0.37266 0.35106 0.65319 0.20673 0.62083 0.39976 0.49995 0.61237 0.12207 0.73728 0.20377 0.18204 0.48058 0.55654 0.63882 0.75948 0.69853 0.69617 0.36073 0.14240 0.42620 0.53601 0.99710 0.37387 0.50305 0.27660 0.77772 0.59910 0.33518 0.12398 0.50510 0.44141 0.62325 0.34660 0.36668 0.54259 0.82641 0.30388 0.62527 0.32578 0.56812 0.81236 0.57249 0.55328 0.44323 0.47118 0.78169 0.57171 0.61229 0.30209 0.43339 0.33588 0.28728 0.16789 0.44649 0.55337 0.21591 0.57889 0.31863 0.39047 0.34375 0.18879 0.19733 0.55145 0.86693 0.52564 0.35898 0.68151 0.37897 0.60868 0.57092 0.72530 0.34375 0.47883 APPENDIX B: TERMINAL OUTPUT FROM SMOOTH ernie> seg smooth Default seed is 123457; do you wish to change it? n DATA FILE (D) OR GENERATE DATA (G)? d ENTER NAME OF DATA FILE. exp25.trans HOW MANY OBSERVATIONS ARE IN exp25.trans ? 100 MAXFN IS BY DEFAULT 30.00000 DO YOU WISH TO CHANGE IT? n EPS IS BY DEFAULT 0.00010 DO YOU WISH TO CHANGE IT? n NSIG IS BY DEFAULT 6.00000 DO YOU WISH TO CHANGE IT? n RHO IS BY DEFAULT 3.00000 DO YOU WISH TO CHANGE IT? n DO YOU WANT TO INPUT THE NUMBER OF INTERVALS OR HAVE AIC OR CAIC COMPUTE IT? (I,A,C) a INTERVAL LIKELIHOOD CRITERION 3 16.84432252 14.84432252 4 21.95344030 18.95344030 5 28.05494509 24.05494509 6 28.01191614 23.01191614 7 28.13292080 22.13292080 8 31.00902772 24.00902772 9 30.71464017 22.71464017 10 30.44856563 21.44856563 11 30.41067909 20.41067909 12 31.08664152 20.08664152 13 33.90164548 21.90164548 14 35.14252143 22.14252143 15 33.65939917 19.65939917 16 33.22997164 18.22997164 17 35.20643562 19.20643562 18 34.00957593 17.00957593 19 32.61787159 14.61787159 20 33.52288774 14.52288774 21 34.67427683 14.67427683 22 36.37304565 15.37304565 23 38.71996709 16.71996709 24 39.24354221 16.24354221 25 38.51366639 14.51366639 26 38.35359312 13.35359312 27 38.82539722 12.82539722 28 38.48174400 11.48174400 29 39.22223041 11.22223041 30 37.83442583 8.83442583 DO YOU WANT TO CHECK FOR MERGING? (Y/N) y CHECKING INTERVALS FOR POSSIBLE MERGING INTERVALS COUNTS DIFFERENCE 1 0.11257 0.31125 0.40285 20 20 -1.92600 2 0.31125 0.40285 0.55333 20 20 -0.21972 3 0.40285 0.55333 0.63204 20 20 -1.06316 4 0.55333 0.63204 0.99710 20 20 -9.76791 NO FURTHER IMPROVEMENT BY MERGING INTERVALS I KNOT COUNT DENSITY CUMULATIVE 1 0.11257 0.00664 2 0.31125 20.00000 1.00662 0.19601 3 0.40285 20.00000 2.18353 0.39535 4 0.55333 20.00000 1.32912 0.59468 5 0.63204 20.00000 2.54065 0.79402 6 0.99710 20.00000 0.54786 0.99336 COEFFICIENTS OF THE SPLINE ARE: 0 0.01512 1 0.05123 2 -0.70768 3 1.62272 4 -0.85474 5 2.33537 6 -1.59243 7 0.13043 RUN NUMBER 1 CURVATURE IS : 123.1605069137 RUN NUMBER IS 1 VALUE AT RIGHT END IS 1.00000 NEW WEIGHTS ARE 0.86 0.86 0.86 0.87 0.87 0.87 0.87 0.88 0.88 0.88 0.88 0.88 0.89 0.89 0.89 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.93 0.93 0.93 0.94 0.94 0.94 0.95 0.95 0.95 0.95 0.95 0.96 0.96 0.96 0.97 0.98 0.98 0.99 0.99 1.00 1.00 1.00 1.00 1.01 1.01 1.01 1.02 1.02 1.02 1.02 1.03 1.03 1.03 1.04 1.04 1.04 1.04 1.04 1.05 1.05 1.05 1.05 1.07 1.07 1.08 1.08 1.09 1.10 1.12 1.12 1.12 1.15 1.16 1.19 1.21 1.22 1.26 1.26 1.28 1.35 1.59 1.64 COEFFICIENTS OF THE SPLINE ARE: 0 0.01666 1 0.04394 2 -0.65359 3 1.54097 4 -0.79860 5 2.24817 6 -1.51693 7 0.11762 RUN NUMBER 2 CURVATURE IS : 108.5770899003 RUN NUMBER IS 2 NEW WEIGHTS ARE 0.86 0.86 0.86 0.87 0.87 0.87 0.87 0.88 0.88 0.88 0.88 0.88 0.89 0.89 0.89 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.93 0.93 0.93 0.94 0.94 0.94 0.95 0.95 0.95 0.95 0.95 0.96 0.96 0.96 0.97 0.98 0.98 0.99 0.99 1.00 1.00 1.00 1.00 1.01 1.01 1.01 1.02 1.02 1.02 1.02 1.03 1.03 1.03 1.04 1.04 1.04 1.04 1.04 1.05 1.05 1.05 1.05 1.07 1.07 1.08 1.08 1.09 1.10 1.12 1.12 1.12 1.15 1.16 1.19 1.21 1.22 1.26 1.26 1.28 1.35 1.59 1.64 INPUT THE NUMBER OF BOOTSTRAP SAMPLES YOU WANT. MAXIMUM = 215 100 INPUT THE SIZE OF EACH BOOTSTRAP SAMPLE. MAXIMUM = 55 1 INPUT NAME OF OUTPUT FILE. exp25_trn.out **** STOP 1 1 APPENDIX C: FULL OUTPUT FROM SMOOTH CUMGEN SCRATCH FILE INTERVAL LIKELIHOOD CRITERION 3 16.84432252 14.84432252 4 21.95344030 18.95344030 5 28.05494509 24.05494509 6 28.01191614 23.01191614 7 28.13292080 22.13292080 8 31.00902772 24.00902772 9 30.71464017 22.71464017 10 30.44856563 21.44856563 11 30.41067909 20.41067909 12 31.08664152 20.08664152 13 33.90164548 21.90164548 14 35.14252143 22.14252143 15 33.65939917 19.65939917 16 33.22997164 18.22997164 17 35.20643562 19.20643562 18 34.00957593 17.00957593 19 32.61787159 14.61787159 20 33.52288774 14.52288774 21 34.67427683 14.67427683 22 36.37304565 15.37304565 23 38.71996709 16.71996709 24 39.24354221 16.24354221 25 38.51366639 14.51366639 26 38.35359312 13.35359312 27 38.82539722 12.82539722 28 38.48174400 11.48174400 29 39.22223041 11.22223041 30 37.83442583 8.83442583 CHECKING INTERVALS FOR POSSIBLE MERGING INTERVALS COUNTS DIFFERENCE 1 0.11257 0.31125 0.40285 20 20 -1.92600 2 0.31125 0.40285 0.55333 20 20 -0.21972 3 0.40285 0.55333 0.63204 20 20 -1.06316 4 0.55333 0.63204 0.99710 20 20 -9.76791 NO FURTHER IMPROVEMENT BY MERGING INTERVALS I KNOT COUNT DENSITY CUMULATIVE 1 0.11257 0.00664 2 0.31125 20.00000 1.00662 0.19601 3 0.40285 20.00000 2.18353 0.39535 4 0.55333 20.00000 1.32912 0.59468 5 0.63204 20.00000 2.54065 0.79402 6 0.99710 20.00000 0.54786 0.99336 1FIXSPL SCRATCH FILE RUN NUMBER IS 1 VALUE AT RIGHT END IS 1.00000 NEW WEIGHTS ARE 0.86 0.86 0.86 0.87 0.87 0.87 0.87 0.88 0.88 0.88 0.88 0.88 0.89 0.89 0.89 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.93 0.93 0.93 0.94 0.94 0.94 0.95 0.95 0.95 0.95 0.95 0.96 0.96 0.96 0.97 0.98 0.98 0.99 0.99 1.00 1.00 1.00 1.00 1.01 1.01 1.01 1.02 1.02 1.02 1.02 1.03 1.03 1.03 1.04 1.04 1.04 1.04 1.04 1.05 1.05 1.05 1.05 1.07 1.07 1.08 1.08 1.09 1.10 1.12 1.12 1.12 1.15 1.16 1.19 1.21 1.22 1.26 1.26 1.28 1.35 1.59 1.64 RUN NUMBER IS 2 NEW WEIGHTS ARE 0.86 0.86 0.86 0.87 0.87 0.87 0.87 0.88 0.88 0.88 0.88 0.88 0.89 0.89 0.89 0.90 0.90 0.90 0.90 0.90 0.90 0.90 0.91 0.91 0.91 0.91 0.91 0.91 0.91 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.92 0.93 0.93 0.93 0.94 0.94 0.94 0.95 0.95 0.95 0.95 0.95 0.96 0.96 0.96 0.97 0.98 0.98 0.99 0.99 1.00 1.00 1.00 1.00 1.01 1.01 1.01 1.02 1.02 1.02 1.02 1.03 1.03 1.03 1.04 1.04 1.04 1.04 1.04 1.05 1.05 1.05 1.05 1.07 1.07 1.08 1.08 1.09 1.10 1.12 1.12 1.12 1.15 1.16 1.19 1.21 1.22 1.26 1.26 1.28 1.35 1.59 1.64 1PLOT SCRATCH FILE COEFFICIENTS OF THE SPLINE ARE: 0 0.01512 1 0.05123 2 -0.70768 3 1.62272 4 -0.85474 5 2.33537 6 -1.59243 7 0.13043 RUN NUMBER 1 CURVATURE IS : 123.1605069137 1 SPLINE ESTIMATE OF CUMULATIVE 0.10E+01 +*********+*********+*********+*********+*********C*********+ * CC * * CCCC * * CCCCCCCC * * CC * 0.90E+00 + C + * C * * C * * C * * C * 0.80E+00 + + * C * * * * C * * C * 0.70E+00 + + * C * * C * * C * * C * 0.60E+00 + C + * C * * * * C * * C * 0.50E+00 + C + * * * C * * C * * C * 0.40E+00 + C + * * * C * * C * * C * 0.30E+00 + C + * C * * * * C * * * 0.20E+00 + C + * C * * C * * * * CC * 0.10E+00 + C + * CC * * C * * CC * * C * 0.00E+00 +*********+*********+*********+*********+*********+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 1 SPLINE ESTIMATE OF DENSITY 0.40E+01 +*********+*********+*********+*********+*********+*********+ * * * * * * * * 0.36E+01 + + * * * * * * * * 0.32E+01 + + * * * * * * * * 0.28E+01 + + * * * * * * * * 0.24E+01 + + * D * * D * * D * * D * 0.20E+01 + D + * D D * * * * D DDDD D * * D DD D * 0.16E+01 + D D D + * DD DDD * * D * * D * * * 0.12E+01 + D + * D * * D * * D D * * D * 0.80E+00 + D D D + * D * * DD D D D * * DD * * D D * 0.40E+00 + D + * D * * D D * * D D * * D D * 0.00E+00 +*********+*********+*********+*********+*DDD*****+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 1 SPLINE ESTIMATE OF SECOND DERIVATIVE 0.25E+02 +*********+*********+*********+*********+*********+*********+ * 2 * * 2 * * * * 2 * 0.20E+02 + + * * * 2 2 * * * * 2 * 0.15E+02 + 2 + * * * 2 2 * * 2 2 * * 2 2 * 0.10E+02 + 2 + * 2 2 * * * * 2 * * 2 2 2 * 0.50E+01 + + * 2 2 2 * * 2 * * 2 * * 2 2 2 * 0.00E+00 +---------------------22----2-------------2-----------------+ * 2 2 * * 2 2 2 * * 2 22 2 * * 2 * -0.50E+01 + 2 2 + * 2 * * 2 * * 2 2 * * 2 * -0.10E+02 + 2 + * 2 * * * * 2 * * 2 * -0.15E+02 + 2 + * 2 * * 2 * * * * 2 * -0.20E+02 + 2 + * 2 * * 2 * * * * * -0.25E+02 +*********+*********+*********+*********+*********+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 1 COEFFICIENTS OF THE SPLINE ARE: 0 0.01666 1 0.04394 2 -0.65359 3 1.54097 4 -0.79860 5 2.24817 6 -1.51693 7 0.11762 RUN NUMBER 2 CURVATURE IS : 108.5770899003 1 SPLINE ESTIMATE OF CUMULATIVE 0.10E+01 +*********+*********+*********+*********+*********C*********+ * CC * * CCCCCCC * * CCCCC * * CC * 0.90E+00 + C + * C * * C * * C * * C * 0.80E+00 + + * C * * * * C * * C * 0.70E+00 + + * C * * C * * C * * C * 0.60E+00 + C + * C * * * * C * * C * 0.50E+00 + C + * * * C * * C * * C * 0.40E+00 + + * C * * C * * C * * C * 0.30E+00 + C + * C * * * * C * * * 0.20E+00 + C + * C * * C * * C * * C * 0.10E+00 + C + * CC * * C * * CC * * C * 0.00E+00 +*********+*********+*********+*********+*********+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 1 SPLINE ESTIMATE OF DENSITY 0.40E+01 +*********+*********+*********+*********+*********+*********+ * * * * * * * * 0.36E+01 + + * * * * * * * * 0.32E+01 + + * * * * * * * * 0.28E+01 + + * * * * * * * * 0.24E+01 + + * * * D * * D D * * * 0.20E+01 + D D + * D * * D * * D DDDD * * D DD D D * 0.16E+01 + D D D D + * DD DDD * * D * * D * * * 0.12E+01 + D D + * * * * * D D D * * * 0.80E+00 + D D D + * D * * D D D D * * DDDD * * D D * 0.40E+00 + + * D D * * D D * * D D * * DD DD * 0.00E+00 +*********+*********+*********+*********+*DD******+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 1 SPLINE ESTIMATE OF SECOND DERIVATIVE 0.25E+02 +*********+*********+*********+*********+*********+*********+ * * * * * 2 * * 2 * 0.20E+02 + + * 2 * * * * 2 * * 2 * 0.15E+02 + 2 + * * * 2 2 * * 2 * * 2 2 * 0.10E+02 + 2 + * 2 2 2 * * 2 * * 2 * * 2 * 0.50E+01 + 2 2 2 + * * * 2 2 2 2 * * 2 2 * * 2 2 2 * 0.00E+00 +----------------------2----2-------------------------------+ * 2 2 2 * * 2 2 * * 2 2 * * 2 2 2 * -0.50E+01 + 2 2 + * 2 2 * * * * 2 2 * * 2 * -0.10E+02 + + * 2 * * 2 * * * * 2 2 * -0.15E+02 + 2 + * 2 * * * * 2 * * 2 * -0.20E+02 + 2 + * 2 * * * * * * * -0.25E+02 +*********+*********+*********+*********+*********+*********+ 0.00E+00 0.40E+00 0.80E+00 0.12E+01 APPENDIX D: THE FORTRAN PROGRAM SMOOTH PROGRAM SMOOTH cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c THIS PROG GETS DATA AND COMPARES THE CALCULATION OF A c c STATISTIC BY BOOTSTRAPPING SAMPLES AND BY A SMOOTHED VERSION c c WHICH APPROXIMATES THE CUMULATIVE USING MONOTONE SPLINES PER c c RAMSEY (STATISTICAL SCIENCE, NOV 1988.) THE MONOTONE SPLINING c c GAYLORD AND RAMIREZ DOES NOT ASSUME POSITIVE COEFFICIENTS OF THE c c B-SPLINE. THE BASIS IS PER RAMSEY, HOWEVER. THIS PROGRAM HAS c c BEEN DESIGNED AND WRITTEN BY GAYLORD AND RAMIREZ, SPRING 1989. c c THE NUMBER OF INTERVALS (KNOTS) USED IS DETERMINED BY USING AIC c c OR INPUT BY USER. TO TRY TO ACHEIVE MONOTONICITY, WEIGHTING IS c c USED ON THE CUMULATIVE. THE SPOTS WHERE MONOTONICITY IS A NOT c c OCCURING GET WEIGHTED MORE. THESE WEIGHTS ARE EXPONENTIALLY c c SMOOTHED IN THIS NEIGHBORHOOD. c c v.6.29.89 c cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER LDX, LDK, LDS, LDBT, LDA PARAMETER ( LDX = 205, + LDK = 35, + LDA = 37, + LDS = 55, + LDBT = 215 ) c THE DIM OF WORK IS MAX(LDA,LDS) DOUBLE PRECISION X(LDX), KNOTS(LDK), F(LDK), + SAM(LDS), SMBT(LDBT), BOOT(LDBT), + EPS, STATCP, RHO, + WORK(LDS), CUMF(LDX), WEIGHT(LDX), + VAL, A(0:LDA), LEFT, + SVAL(LDX, LDA),SS((LDA+1)*(LDA+2)/2), + RIGHT, TOP DOUBLE PRECISION DSEED, START INTEGER NDATA, NKNOT, SAMPNO, + MAXSAM, SMPSZE, IJOB, + IER, INTERV, NSIG, + MAXFN, IWK(2*LDA), RUN, + IERLEQ, CNT(LDK) LOGICAL FIX EXTERNAL GETDAT, GENBT, GENSMD, + CUMGEN, STATCP, SDGEN, + OUTPUT, GETSIZ, ENDPTS, + EPSSET, BSVALS, FIXSPL, + VTPROF, VMULFM, LEQ1S, + FCOMP, PLOT, CUMFVL, + SETWT INTRINSIC DBLE, REAL, SQRT COMMON /SEED/ DSEED + /COEF/ A + /WORK/ WORK + /DATA/ NDATA, X, WEIGHT + /KNOTS/ NKNOT, KNOTS, LEFT, RIGHT, TOP + /INDEX/ INTERV + /VALUE/ VAL, + /COUNTS/ CNT DATA NSIG /6/ + MAXFN /30/ + EPS /1.0D-4/ + RHO /3.0D0/ c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c VARIABLE NAMES: c c X ORIGINAL DATA c c CUMF EMPIRICAL CUMULATIVE DISTRIBUTION AT X c c NDATA NUMBER OF ORIGINAL DATA c c KNOTS ABSISSAE OF CUMULATIVE c c F ORDINATES OF CUMULATIVE c c NKNOT NUMBER OF POINTS IN KNOTS AND F c c A COEFFICIENTS OF B-SPLINE. c c WORK NECESSARY WORK SPACE FOR IMSL c c SAM SAMPLE OF BOOTSTRAPPED DATA; OVERWRITTEN OFTEN c c SMPSZE SIZE OF BOOTSTRAPPED SAMPLE c c SMBT SAMPLE OF SMOOTHED STATISTIC c c BOOT SAMPLE OF USUAL BOOTSTAPPED STATISTIC c c SAMPNO COUNTER FOR SMBT AND BOOT c c MAXSAM MAXIMUM NUMBER OF SAMPLES; UPPER LIMIT OF SAMPNO c c LDX,LDF, c c LDS,LDBT, c c LDA DIMENSIONS OF RESPECTIVE ARRAYS c c c c ROUTINE/FUNCTION NAMES: c c SDGEN ROUTINE TO START SEED c c GETDAT ROUTINE TO GET DATA c c CUMGEN ROUTINE TO GENERATE CUMULATIVE c c SPLSET ROUTINE TO DETERMINE A, THE COEFFICIENTS OF S c c GENSMD ROUTINE TO GENERATE SMOOTHED BOOTSTRAP SAMPLE c c GENBT ROUTINE TO GENERATE USUAL BOOTSTRAP SAMPLE c c STATCP FUNCTION TO CALCULATE SAMPLE STATISTIC c c OUTPUT OUTPUT ROUTINE c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c RUN = 0 CALL SDGEN START = DSEED CALL GETDAT (X, NDATA) CALL EPSSET (MAXFN, EPS, NSIG, RHO) OPEN (UNIT = 7, FILE = 'JUNK.CUMGEN') CALL CUMGEN (X, KNOTS, F, NDATA, NKNOT) CLOSE (UNIT = 7) CALL SETWT (WEIGHT, NDATA) OPEN (UNIT = 8, FILE = 'JUNK.PLOT') WRITE (8, 81) 'PLOT SCRATCH FILE' OPEN (UNIT = 9, FILE = 'JUNK.FIXSPL') WRITE (9, 81) 'FIXSPL SCRATCH FILE' 50 CONTINUE VAL = 0.0D0 RUN = RUN + 1 CALL BSVALS (X, KNOTS, WEIGHT, NDATA, NKNOT, SVAL) CALL CUMFVL (CUMF, WEIGHT, NDATA) CALL VTPROF (SVAL, NDATA, NKNOT+2, LDX, SS) CALL VMULFM (SVAL, CUMF, NDATA, NKNOT+2, 1, LDX, LDX, A, LDA, IER) IJOB = 0 IERLEQ = 0 CALL LEQ1S (SS, NKNOT+2, A, 1, LDA, IJOB, IWK, WORK, IERLEQ) CALL PLOT(X, KNOTS, NDATA, NKNOT, RUN) IF (IERLEQ .EQ. 0) THEN CALL FIXSPL (FIX, RUN, RHO) IF (FIX) GOTO 50 ENDIF CLOSE (UNIT = 8) CLOSE (UNIT = 9) CALL ENDPTS CALL FCOMP (KNOTS, F, NKNOT) CALL GETSIZ (MAXSAM, SMPSZE, LDS, LDBT) DO 100 SAMPNO = 1,MAXSAM CALL GENSMD (F, NSIG, MAXFN, EPS, SAM, SMPSZE, WORK) SMBT (SAMPNO) = STATCP (SAM, SMPSZE, LDS) CALL GENBT (X, NDATA, SAM, SMPSZE) BOOT (SAMPNO) = STATCP (SAM, SMPSZE, LDS) 100 CONTINUE CALL OUTPUT (X, F, NDATA, SMBT, BOOT, SMPSZE, + MAXSAM, START, A) 81 FORMAT (/1X, A, /) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE SDGEN c c THIS ROUTINE WILL GENERATE THE SEED OF THE PROGRAM c THERE IS AN OPTION OF INPUTTING THE SEED OR LETTING IT BE CHOSEN c DOUBLE PRECISION DSEED INTRINSIC DBLE CHARACTER*1 LETTER COMMON /SEED/ DSEED DSEED = DBLE(123457) WRITE (*,*) ' Default seed is 123457; do you wish to change it?' READ (*,'(A)') LETTER IF ((LETTER .EQ. 'Y') .OR. (LETTER .EQ. 'y')) THEN WRITE (*,*) ' INPUT SEED' READ (*,*) DSEED ENDIF RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE GETDAT (X, N) INTEGER N, LDX, I PARAMETER (LDX = 205) DOUBLE PRECISION X(LDX), LAMBDA, DSEED CHARACTER*1 ANS CHARACTER*80 FILENM LOGICAL EX COMMON /SEED/ DSEED EXTERNAL GGEXN 10 WRITE (*,*) ' DATA FILE (D) OR GENERATE DATA (G)?' READ (*,199) ANS IF ((ANS .EQ. 'G') .OR. (ANS .EQ.'g')) THEN WRITE (*,*) ' DATA WILL BE GENERATED BY AN EXPONENTIAL DIST.' WRITE (*,*) ' INPUT LAMBDA, MEAN OF DIST.' READ (*,*) LAMBDA WRITE (*,*) ' HOW MANY DEVIATES? (MAX =',LDX,')' READ (*,*) N CALL GGEXN (DSEED, LAMBDA, N, X) ELSE WRITE (*,*) ' ENTER NAME OF DATA FILE.' READ (*,199) FILENM INQUIRE (FILE = FILENM, EXIST = EX) IF (.NOT. EX) THEN WRITE (*,*) ' FILE ',FILENM,' DOES NOT EXIST!' GOTO 10 ENDIF WRITE (*,*) ' HOW MANY OBSERVATIONS ARE IN ',FILENM,'?' READ (*,*) N OPEN (FILE = FILENM, UNIT = 5, ERR = 9999) DO 100 I=1,N READ (UNIT=5,FMT=*,ERR=9999) X(I) 100 CONTINUE CLOSE (UNIT = 5) ENDIF RETURN 199 FORMAT (A) 9999 WRITE (*,*) ' ERROR IN OPENING OR READING FROM ',FILENM GOTO 10 END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE EPSSET (MAXFN, EPS, NSIG, RHO) DOUBLE PRECISION EPS, RHO REAL VAL(4) INTEGER NSIG, I, MAXFN INTRINSIC REAL, INT, DBLE CHARACTER*6 LABEL(4) CHARACTER*1 ANS DATA LABEL / ' MAXFN', ' EPS ', ' NSIG', ' RHO'/ VAL(1) = REAL(MAXFN) VAL(2) = REAL(EPS) VAL(3) = REAL(NSIG) VAL(4) = REAL(RHO) DO 15 I=1,4 WRITE (*, 80) LABEL(I),' IS BY DEFAULT ',VAL(I) WRITE (*, *) ' DO YOU WISH TO CHANGE IT?' READ (*, 100) ANS IF ((ANS .EQ.'Y') .OR. (ANS .EQ. 'y')) THEN WRITE (*, *) ' INPUT ',LABEL(I) READ (*, *) VAL(I) ENDIF 15 CONTINUE MAXFN = INT(VAL(1)) EPS = DBLE(VAL(2)) NSIG = INT(VAL(3)) RHO = DBLE(VAL(4)) RETURN 80 FORMAT (1X, A6, A15, F10.5) 100 FORMAT (A1) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE CUMGEN (X, X1, F, N, M) ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc c THIS SUBROUTINE GROUPS DATA INTO CLASS INTERVALS BY USING A MODEL c c SELECTION HISTOGRAM AND CREATING AN EMPIRICAL CUMULATIVE FUNCTION c c USING THESE CLASS INTERVALS. THE INTERVALS ARE EVENLY WEIGHTED. c c APRIL 23, 1989 c ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc INTEGER LDK, LX PARAMETER ( LDK = 35, + LX = 205 ) DOUBLE PRECISION X(LX), X1(LDK), F(LDK), AIC, TMPX(LDK), + ZERO, LOGN, SET, CAIC, DENS(LDK), + LIKEHD,DCRITV, CRITVL, TMPF(LDK), + DD, NN, ONE, DN, + TMPDEN(LDK), THIRD INTEGER COUNT, REM, STOPIT, CRIT, START, + MM, N, M, I, J, + CNT(LDK),TMPCNT(LDK) LOGICAL MERGE CHARACTER*1 ANS COMMON /COUNTS/ CNT EXTERNAL VSRTAD, DETSTP INTRINSIC MOD, LOG, DBLE DATA ONE /1.0D0/, + ZERO /0.0D0/ c SET = -1.0 D38 DN = DBLE(N) LOGN = LOG(DN) THIRD = ONE / (3.0D0) cALL VSRTAD (X, N) c c CRIT = 1 => USER INPUT NUMBER OF INTERVALS c 2 => USE AIC TO DETERMINE NUMBER OF INTERVALS c 3 => USE CAIC TO DETERMINE NUMBER OF INTERVALS c CALL DETSTP (CRIT, START, STOPIT, N, LDK) WRITE (7,69) 'CUMGEN SCRATCH FILE', + 'INTERVAL','LIKELIHOOD','CRITERION' WRITE (*,71) 'INTERVAL','LIKELIHOOD','CRITERION' DO 200 I = START, STOPIT LIKEHD = ZERO REM = MOD (N,I) DO 10 J=1,REM TMPCNT(J) = ONE + N/I 10 CONTINUE DO 20 J=REM+1,I TMPCNT(J) = N/I 20 CONTINUE COUNT = 0 TMPX (1) = X (1) TMPF (1) = (ONE - THIRD) / (DN + THIRD) c ADD EXTRA VALUE AT END FOR AVERAGING X(N+1) = X(N) DO 30 J=2,I+1 COUNT = COUNT + TMPCNT(J-1) TMPF (J) = (COUNT - THIRD) / (DN + THIRD) TMPX (J) = (X(COUNT) + X(COUNT+1))/2.0 TMPDEN (J-1) = TMPCNT(J-1) / (TMPX(J)-TMPX(J-1)) LIKEHD = LIKEHD + TMPCNT(J-1) * LOG(TMPDEN(J-1)) 30 CONTINUE LIKEHD = LIKEHD -DN*LOGN IF (CRIT .EQ. 2) THEN AIC = LIKEHD - DBLE(I-1) CRITVL = AIC ELSE IF (CRIT .EQ. 3) THEN CAIC = LIKEHD - LOGN * 0.5D0 * DBLE(I-1) CRITVL = CAIC ENDIF WRITE (7,72) I,LIKEHD,CRITVL WRITE (*,72) I,LIKEHD,CRITVL IF (CRITVL .GT. SET) THEN SET = CRITVL M = I + 1 DO 40 J=1,M F(J) = TMPF(J) X1(J) = TMPX(J) 40 CONTINUE DO 45 J=1,M-1 CNT(J) = TMPCNT(J) DENS(J) = TMPDEN(J) 45 CONTINUE ENDIF 200 CONTINUE CRITVL = SET MERGE = .FALSE. IF (CRIT .NE. 1) THEN WRITE (*,*) ' DO YOU WANT TO CHECK FOR MERGING? (Y/N)' READ (*,'(A)') ANS MERGE = ((ANS .EQ. 'Y') .OR. (ANS .EQ. 'y')) ENDIF IF ((CRIT .NE. 1) .AND. MERGE) THEN c CHECK FOR POSSIBLE MERGING OF INTERVALS 210 CONTINUE WRITE (7,73) 'CHECKING INTERVALS FOR POSSIBLE MERGING', + 'INTERVALS','COUNTS','DIFFERENCE' WRITE (*,73) 'CHECKING INTERVALS FOR POSSIBLE MERGING', + 'INTERVALS','COUNTS','DIFFERENCE' SET = -1.0 D38 DO 50, J = 2, M-1 NN = CNT(J-1)+CNT(J) DD = NN / (X1(J+1)-X1(J-1)) DCRITV = NN * LOG(DD) - - CNT(J-1) * LOG(DENS(J-1)) - - CNT(J) * LOG(DENS(J)) IF (CRIT .EQ. 2) THEN DCRITV = DCRITV - (1.0D0 - 2.0D0) ELSE DCRITV = DCRITV - 0.5D0 * * (LOG(DBLE(M-3)) - LOG(DBLE(M-2))) ENDIF WRITE (7,74) J-1,X1(J-1),X1(J),X1(J+1),CNT(J-1),CNT(J), + DCRITV WRITE (*,74) J-1,X1(J-1),X1(J),X1(J+1),CNT(J-1),CNT(J), + DCRITV IF (DCRITV .GT. SET) THEN SET = DCRITV MM = J ENDIF 50 CONTINUE IF (SET .GT. ZERO) THEN CRITVL = CRITVL + SET M = M - 1 WRITE (7,75) 'IMPROVEMENT POSSIBLE', + 'MERGING INTERVALS ',MM-1,' AND ',MM, + 'NEW CRITERION VALUE IS ',CRITVL WRITE (*,75) 'IMPROVEMENT POSSIBLE', + 'MERGING INTERVALS ',MM-1,' AND ',MM, + 'NEW CRITERION VALUE IS ',CRITVL CNT(MM-1) = CNT(MM-1) + CNT(MM) DENS(MM-1)= CNT(MM-1) / (X1(MM+1)-X1(MM-1)) DO 60, J = MM,M F(J) = F(J+1) X1(J) = X1(J+1) 60 CONTINUE DO 70, J = MM,M-1 CNT(J) = CNT(J+1) DENS(J) = DENS(J+1) 70 CONTINUE IF (M .LE. 4) THEN GOTO 273 ELSE GOTO 210 ENDIF ELSE WRITE (7,76) 'NO FURTHER IMPROVEMENT BY MERGING INTERVALS' WRITE (*,76) 'NO FURTHER IMPROVEMENT BY MERGING INTERVALS' ENDIF ENDIF WRITE (7, 81) 'I', 'KNOT', 'COUNT', 'DENSITY', 'CUMULATIVE' WRITE (7, 82) 1, X1(1), F(1) WRITE (*, 81) 'I', 'KNOT', 'COUNT', 'DENSITY', 'CUMULATIVE' WRITE (*, 82) 1, X1(1), F(1) DO 250 I = 2, M WRITE (7, 83) I, X1(I), CNT(I-1), DENS(I-1)/DN, F(I) WRITE (*, 83) I, X1(I), CNT(I-1), DENS(I-1)/DN, F(I) 250 CONTINUE 273 CONTINUE RETURN 69 FORMAT (1X, A // 1X, A9, A20, A20 /) 71 FORMAT (1X, // 1X, A9, A20, A20 /) 72 FORMAT (1X, I9, F20.8, F20.8) 73 FORMAT (/1X, A / 1X, A35, A10, A15 //) 74 FORMAT (1X, I5, F10.5, F10.5, F10.5, I5, I5, F15.5) 75 FORMAT (1X, A / 1X, A, I9, A, I9 / 1X, A, F20.8 /) 76 FORMAT (1X, A) 81 FORMAT (/ 1X, A5, A12, A12, A12, A12 /) 82 FORMAT ( 1X, I5, F12.5, 24X, F12.5) 83 FORMAT ( 1X, I5, 4(F12.5)) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE DETSTP (CRIT, START, STOP, N, LDK) c APRIL 23, 1989 INTEGER CRIT, START, STOP, N, LDK CHARACTER*1 LETTER INTRINSIC SQRT, REAL, MIN, INT, MAX 100 CONTINUE WRITE (*,*) ' DO YOU WANT TO INPUT THE NUMBER OF INTERVALS' WRITE (*,*) ' OR HAVE AIC OR CAIC COMPUTE IT? (I,A,C)' READ (*,1000) LETTER IF ((LETTER .EQ. 'I') .OR. (LETTER .EQ. 'i')) THEN CRIT = 1 WRITE (*,*) ' INPUT NUMBER OF INTERVALS (AT LEAST 4).' READ (*,*) START STOP = START ELSE IF ((LETTER .EQ. 'a') .OR. (LETTER .EQ. 'A')) THEN CRIT = 2 STOP = MIN (2*INT(SQRT(REAL(N))) + 10, N/2) START = 3 ELSE IF ((LETTER .EQ. 'C') .OR. (LETTER .EQ. 'c')) THEN CRIT = 3 STOP = MIN (2*INT(SQRT(REAL(N))) + 5, N/2) START = 3 ELSE GOTO 100 ENDIF STOP = MIN (STOP, LDK-1) START = MAX (3, START) RETURN 1000 FORMAT (A) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE SETWT (WEIGHT, NDATA) INTEGER LDX, I, NDATA PARAMETER (LDX = 205) DOUBLE PRECISION ONE, WEIGHT(LDX) INTRINSIC DBLE ONE = 1.0D0 DO 100, I = 1, NDATA WEIGHT(I) = ONE 100 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE BSVALS (X, T, WEIGHT, NX, NK, SVAL) INTEGER LDX, NX, NK, LDK, LDA, + I, K, INT PARAMETER (LDX = 205, LDK = 35, LDA = 37) DOUBLE PRECISION X(LDX), T(LDK), SVAL(LDX,0:LDA), BS, + WTI, WEIGHT(LDX) EXTERNAL BS, FINDIT INTRINSIC SQRT COMMON /INDEX/ INT DO 100 I = 1, NX WTI = SQRT (WEIGHT(I)) CALL FINDIT (X(I), T, NK, INT, LDK) DO 150 K = 0, NK+1 SVAL(I,K) = BS (X(I),K) * WTI 150 CONTINUE 100 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE CUMFVL (CUMF, WEIGHT, NDATA) INTEGER LDX, NDATA, I PARAMETER (LDX = 205) DOUBLE PRECISION CUMF(LDX), WEIGHT(LDX), THIRD INTRINSIC DBLE, SQRT THIRD = 1.0D0 / 3.0D0 DO 50 I = 1, NDATA CUMF(I) = DBLE(I-THIRD) / DBLE(NDATA+THIRD) * SQRT(WEIGHT(I)) 50 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE PLOT (X, T, NDATA, NABSI, RUN) INTEGER LDX, NDATA, NINT, I, NXLABL, NYLABL, + NABSI, IOPT, NTITLE, IXLABL, IYLABL, IER, + LDK, RUN, LDA PARAMETER (LDX = 205, + LDK = 35, + LDA = 36, + NINT = 59) DOUBLE PRECISION XVAL, GAP, CURVTR, X(LDX), T(LDK), + A(0:LDA) REAL RANGE(4),RLZERO, PTWORK(0:NINT,4), S, DS, DDS CHARACTER ICHAR*1, ITITLE*72 EXTERNAL S, DS, DDS, USPLO, UGETIO INTRINSIC REAL COMMON /COEF/ A RLZERO = 0.0E0 NTITLE = 72 WRITE (8,68) '1 COEFFICIENTS OF THE SPLINE ARE:' WRITE (*,*) ' COEFFICIENTS OF THE SPLINE ARE:' DO 50 I = 0, NABSI+1 WRITE (8,75) I, A(I) WRITE (*,75) I, A(I) 50 CONTINUE GAP = (X(NDATA) - X(1)) / NINT CURVTR = RLZERO DO 80, I = 0, NINT XVAL = X(1) + I*GAP PTWORK(I,2) = S(XVAL) PTWORK(I,3) = DS(XVAL) PTWORK(I,4) = DDS(XVAL) PTWORK(I,1) = REAL(XVAL) CURVTR = CURVTR + PTWORK(I,4)**2 80 CONTINUE CURVTR = CURVTR / (NINT + 1.0D0) NXLABL = 0 NYLABL = 0 DO 85 I = 1, 4 RANGE(I) = RLZERO 85 CONTINUE WRITE (*, 81) 'RUN NUMBER ', RUN WRITE (8, 81) 'RUN NUMBER ', RUN WRITE (*,89) 'CURVATURE IS : ', CURVTR WRITE (8,89) 'CURVATURE IS : ', CURVTR CALL UGETIO (3, 0, 8) ITITLE = 'SPLINE ESTIMATE OF CUMULATIVE' ICHAR = 'C' IOPT = 0 CALL USPLO(PTWORK(0,1),PTWORK(0,2),1+NINT,NINT+1,1,1,ITITLE, + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) ITITLE = 'SPLINE ESTIMATE OF DENSITY' ICHAR = 'D' CALL USPLO(PTWORK(0,1),PTWORK(0,3),1+NINT,NINT+1,1,1,ITITLE, + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) ITITLE = 'SPLINE ESTIMATE OF SECOND DERIVATIVE' ICHAR = '2' CALL USPLO(PTWORK(0,1),PTWORK(0,4),1+NINT,NINT+1,1,1,ITITLE, + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) 68 FORMAT (A) 75 FORMAT ( 1X, I5, F10.5) 81 FORMAT (/1X, A, I10, /) 89 FORMAT (/1X, A, 2X, F20.10, /) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE FIXSPL (FIX, RUN, LAMBDA) INTEGER LDK, LDX, RUN PARAMETER (LDK = 35, LDX = 205) DOUBLE PRECISION DDI, DDI1, ONE, ZERO, KNTS(LDK), + EPS, ADDX, TSLOPE, X(LDX), LAMBDA, WEIGHT(LDX), + TMPADD INTEGER NX, CNT(LDK), NKNOTS, I, IADD REAL DS, DMIN, S, DSI, S1, SN, + DDS, D1, D2, D3, RLZERO, DDSI LOGICAL FIX, FIXN, FIX1, FIXM, FIXE EXTERNAL DS, DDS, S, EXPSM INTRINSIC DBLE, SQRT, ABS, MIN, REAL, MAX, EXP COMMON /KNOTS/ NKNOTS, KNTS + /COUNTS/ CNT + /DATA/ NX, X, WEIGHT DATA + ONE /1.0D0/ + ZERO /0.0D0/ + EPS /1.0D-08/ + RLZERO /0.0E0/ WRITE (*, 80) 'RUN NUMBER IS ', RUN WRITE (9, 80) 'RUN NUMBER IS ', RUN FIX = .FALSE. FIXN = .FALSE. FIX1 = .FALSE. FIXM = .FALSE. FIXE = .FALSE. DMIN = RLZERO DO 100 I = 1, NKNOTS-1 DDI = DBLE(DDS (KNTS(I)+EPS)) DDI1 = DBLE(DDS (KNTS(I+1)-EPS)) TSLOPE = (DDI1-DDI) / (KNTS(I+1)-EPS-KNTS(I)-EPS) TMPADD = KNTS(I) - DDI/TSLOPE IF ((TMPADD .GT. KNTS(I)).AND.(TMPADD .LT. KNTS(I+1))) THEN DSI = DS (TMPADD) DDSI= DDS(TMPADD) IF ((DSI .LT. DMIN) .AND. (DDSI .GT. RLZERO)) THEN IADD = I DMIN = DSI ADDX = TMPADD FIX = .TRUE. FIXM = .TRUE. ENDIF ENDIF D1 = DS(KNTS(I)) D2 = DS(KNTS(I+1)) D3 = MIN(D1,D2) IF (D3 .LT. DMIN) THEN IADD = I DMIN = D3 IF (D1 .GT. D3) THEN ADDX = KNTS(I+1) ELSE ADDX = KNTS(I) ENDIF FIX = .TRUE. FIXE = .TRUE. ENDIF 100 CONTINUE SN = S(KNTS(NKNOTS) - EPS) S1 = S(KNTS(1) + EPS) FIXN = (SN.GT.ONE) FIX1 = (S1.LT.ZERO) IF (FIX) THEN WRITE (*, 81) 'MIN IS ', DMIN, ' AT ', ADDX WRITE (9, 81) 'MIN IS ', DMIN, ' AT ', ADDX CALL EXPSM (ADDX, LAMBDA) ENDIF IF (FIXN) THEN FIX = .TRUE. WRITE (*, 81) 'VALUE AT RIGHT END IS ', SN WRITE (9, 81) 'VALUE AT RIGHT END IS ', SN CALL EXPSM (X(NX), LAMBDA) ENDIF IF (FIX1) THEN FIX = .TRUE. WRITE (*, 81) 'VALUE AT LEFT END IS ', S1 WRITE (9, 81) 'VALUE AT LEFT END IS ', S1 CALL EXPSM (X(1), LAMBDA) ENDIF WRITE (*, 81) 'NEW WEIGHTS ARE ' WRITE (*, 83) (WEIGHT(I), I=1, NX) WRITE (9, 81) 'NEW WEIGHTS ARE ' WRITE (9, 83) (WEIGHT(I), I=1, NX) 80 FORMAT (/1X, A, I10) 81 FORMAT (/1X, A, F10.5, A, F10.5) 83 FORMAT ((1X, 10F7.2)) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE EXPSM (PIVOT, LAMBDA) INTEGER I, NX, NKNOTS, LDX, LDF PARAMETER (LDX = 205, LDF = 35) DOUBLE PRECISION PIVOT, KNTS(LDF), X(LDX), WEIGHT(LDX), + LAMBDA, ONE, RANGE, SUM INTRINSIC ABS, EXP, DBLE COMMON /KNOTS/ NKNOTS, KNTS + /DATA/ NX, X, WEIGHT ONE = 1.0D0 SUM = 0.0D0 RANGE = KNTS(NKNOTS) - KNTS(1) DO 100, I = 1, NX WEIGHT(I) = WEIGHT(I) * * (ONE + EXP(-LAMBDA * (ABS (X(I)-PIVOT) / RANGE))) SUM = SUM + WEIGHT(I) 100 CONTINUE c MAKE WEIGHTS ADD TO NX DO 200 I=1,NX WEIGHT(I) = DBLE(NX) * WEIGHT(I) / SUM 200 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE ENDPTS INTEGER NKNOT, LDA, LDF, I PARAMETER (LDA = 37, LDF = 35) DOUBLE PRECISION LEFT, RIGHT, A(0:LDA), ZERO, T(LDF), TOP, ONE, + VAL, SLOPE REAL S, DS EXTERNAL S, DS INTRINSIC DBLE, MAX COMMON /KNOTS/ NKNOT, T, LEFT, RIGHT, TOP + /COEF/ A + /VALUE/ VAL ZERO = 0.0D0 ONE = 1.0D0 VAL = ZERO TOP = ZERO DO 100 I=0,NKNOT+1 TOP = TOP + A(I) 100 CONTINUE c SLOPE = DBLE(DS(T(1))) SLOPE = MAX(SLOPE, 5.0D-3) LEFT = T(1) - S(T(1))/SLOPE c SLOPE = DBLE(DS(T(NKNOT))) SLOPE = MAX(SLOPE, 5.0D-3) RIGHT = T(NKNOT) + (ONE-TOP) / SLOPE c IF (LEFT .GT. T(1)) LEFT = T(1) IF (RIGHT .LT. T(NKNOT)) RIGHT = T(NKNOT) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE FCOMP (X, F, NF) INTEGER NF, LDF, I PARAMETER (LDF = 35) DOUBLE PRECISION F(LDF), X(LDF), VAL REAL S EXTERNAL S COMMON /VALUE/ VAL VAL = 0.0D0 DO 100 I = 1,NF F(I) = S(X(I)) 100 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE GETSIZ (NBOOTS, SZBOOT, LDS, LDBT) INTEGER NBOOTS, SZBOOT, LDS, LDBT INTRINSIC MIN WRITE (*,*) ' INPUT THE NUMBER OF BOOTSTRAP SAMPLES YOU WANT.', + ' MAXIMUM = ', LDBT READ (*,*) NBOOTS NBOOTS = MIN (NBOOTS, LDBT) WRITE (*,*) ' INPUT THE SIZE OF EACH BOOTSTRAP SAMPLE.', + ' MAXIMUM = ', LDS READ (*,*) SZBOOT SZBOOT = MIN (SZBOOT, LDS) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE GENSMD (F, NSIG, ITMAX, EPS, SAMP, K, WORK) c c THIS SUBROUTINE WILL GENERATE A RANDOM SAMPLE FROM AN EMPIRICAL c SMOOTHED DISTRIBUTION. c APRIL 26, 1987 INTEGER LDK, LDS PARAMETER ( LDK = 35, + LDS = 55) DOUBLE PRECISION T(LDK), EPS, F(LDK), T1, T2, + VAL, LEFT, RIGHT, + SAMP(LDS), SMALL DOUBLE PRECISION DSEED REAL S, S1, S2, WORK(LDS) INTEGER K, I, ITMAX, IER, INTERV, + NSIG, NABSI, HOLD INTRINSIC DBLE, ABS EXTERNAL ZFALSE, S, GGUBS, FINDIT COMMON /SEED/ DSEED + /VALUE/ VAL + /KNOTS/ NABSI, T, LEFT, RIGHT HOLD = ITMAX SMALL = 1.0E-16 c GGUBS IS SINGLE PRECISION CALL GGUBS (DSEED, K, WORK) DO 100 I=1,K VAL = DBLE(WORK(I)) SAMP(I) = VAL CALL FINDIT (SAMP(I), F, NABSI, INTERV, LDK) IF (INTERV .EQ. 0) THEN T1 = LEFT T2 = T(1) ELSE IF (INTERV .EQ. NABSI) THEN T1 = T(NABSI) T2 = RIGHT ELSE T1 = T(INTERV) T2 = T(INTERV+1) ENDIF IER = 0 S1 = S(T1) S2 = S(T2) IF (ABS(T1-T2) .LT. SMALL) THEN SAMP(I) = T1 ELSE CALL ZFALSE (S, EPS, NSIG, T1, T2, SAMP(I), ITMAX, IER) ENDIF IF (IER .NE. 0) + SAMP(I) = T1 + (T2-T1) * (SAMP(I)-S1) / (S2-S1) IER = 0 ITMAX = HOLD 100 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE GENBT (X, N, SAMP, M) c c THIS SUBROUTINE WILL GENERATE A BOOTSTRAP SAMPLE FROM X. c INTEGER LDS, LDX PARAMETER ( LDS = 55, + LDX = 205 ) DOUBLE PRECISION X(LDX), SAMP(LDS) INTEGER N, M, I, ISAMP(LDS) DOUBLE PRECISION DSEED COMMON /SEED/ DSEED EXTERNAL GGUD c c GENERATE UNIFORM (DISCRETE) 1,N INDICES. c CALL GGUD (DSEED, N, M, ISAMP) DO 100 I=1,M SAMP(I) = X(ISAMP(I)) 100 CONTINUE RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION STATCP (SAMP, SIZE, LDS) c c THIS ROUTINE WILL CALCULATE THE MEAN OF SAMPLE. c INTEGER SIZE, I, LDS DOUBLE PRECISION SAMP(LDS), STATCP INTRINSIC DBLE STATCP = 0.0D0 DO 100 I=1,SIZE STATCP = STATCP + SAMP(I) 100 CONTINUE STATCP = STATCP/DBLE(SIZE) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE OUTPUT(X, F, NX, SMTHBT, REGBT, SZBT, + NBT, START, A) c INTEGER LDX, LDK, LDA, LDBT PARAMETER ( LDX = 205, + LDK = 35, + LDA = 37, + LDBT = 215 ) DOUBLE PRECISION X(LDX), X1(LDK), F(LDK), SMTHBT(LDBT), + A(0:LDA),START, VAL, REGBT(LDBT), GAP, + RIGHT, XVAL, LEFT REAL S, DS, SI, DSI, DDSI, DDS, + CRTURE INTEGER NX, NX1, NBT, SZBT, INT, I, NINT CHARACTER*20 FILENM CHARACTER*1 ANS LOGICAL FEX COMMON /INDEX/ INT + /VALUE/ VAL + /KNOTS/ NX1, X1, LEFT, RIGHT EXTERNAL S, DS, DDS INTRINSIC DBLE DATA NINT /100/ VAL = 0.0D0 15 WRITE (*,*) ' INPUT NAME OF OUTPUT FILE.' READ (*,1000) FILENM INQUIRE (FILE = FILENM, EXIST = FEX) IF (FEX) THEN WRITE (*,*) ' THAT FILE ALREADY EXISTS;', + ' DO YOU WISH TO OVERWRITE IT? (Y/N)' READ (*,'(A)') ANS IF ((ANS .EQ. 'n') .OR. (ANS .EQ. 'N')) GOTO 15 ENDIF OPEN (FILE = FILENM, UNIT = 6) WRITE (6,*) ' THE ORIGINAL SEED: ',START WRITE (6,*) ' THE ORIGINAL DATA' DO 100 I=1,NX SI = S(X(I)) DSI = DS(X(I)) DDSI = DDS(X(I)) WRITE (6,1500) X(I), SI, DSI, DDSI 100 CONTINUE WRITE (6,*) ' THE CUMULATIVE AND ITS APPROXIMATION AT X1 POINTS' DO 200 I=1,NX1 WRITE (6,2000) X1(I), F(I) 200 CONTINUE WRITE (6,*) ' THE BASIS COEFFICIENTS' DO 205 I=0,NX1+1 WRITE (6,3000) I, A(I) 205 CONTINUE WRITE (6,*) 'THE VALUES FOR PLOTTING' WRITE (6,1400) 'X-VALUE','SPLINE','DENSITY','SECOND DERIV' GAP = RIGHT - LEFT GAP = GAP / NINT CRTURE = 0.0E0 DO 210, I=0, NINT XVAL = LEFT + DBLE(I)*GAP DDSI = DDS(XVAL) WRITE (6,1500) XVAL, S(XVAL), DS(XVAL), DDSI CRTURE = CRTURE + DDSI**2 210 CONTINUE CRTURE = CRTURE * (X(NX) - X(1)) / (NINT + 1) WRITE (6,1100) 'ESTIMATED CURVATURE IS ',CRTURE WRITE (6,*) + ' THE SMOOTHED AND UNSMOOTHED BOOTSTRAPS SAMPLES SIZE ',SZBT DO 300 I=1,NBT WRITE (6,2000) SMTHBT(I),REGBT(I) 300 CONTINUE CLOSE (UNIT = 6) RETURN 1000 FORMAT (A) 1100 FORMAT (/1X, A, F16.8/) 1400 FORMAT (/1X, A16, 3(1X, A16 )) 1500 FORMAT ( 1X, F16.8, 3(1X, F16.8)) 2000 FORMAT ( 1X, F16.8, 2X, F16.8) 3000 FORMAT ( 1X, I5, 3X, F16.8) END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc SUBROUTINE FINDIT (VAL, Y, M, IT, LDY) C C THIS ROUTINE FINDS THE INTERVAL [Y(IT),Y(IT+1)) WHICH CONTAINS C VAL. Y MUST BE STRICTLY INCREASING. THIS VERSION WORKS BY C BISECTION. C INTEGER M, IT, NEW, BOTTOM, TOP, INC, LDY DOUBLE PRECISION Y(LDY), VAL IF (VAL .GT. Y(M)) THEN IT = M RETURN ELSE IF (VAL .LT. Y(1)) THEN IT = 0 RETURN ENDIF NEW = M/2 BOTTOM = 1 TOP = M - 1 INC = 1 10 IF ((VAL .GE. Y(NEW)) .AND. (VAL .LE. Y(NEW+1))) THEN IT = NEW RETURN ELSE IF (VAL .GE. Y(NEW+1)) THEN BOTTOM = NEW + 1 NEW = TOP - (TOP-NEW)/2 INC = INC + 1 ELSE TOP = NEW NEW = BOTTOM + (NEW-BOTTOM)/2 INC = INC + 1 ENDIF GOTO 10 END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION S(X) INTEGER NKNOT, INT, LDK, LDA, I PARAMETER ( LDA = 37, + LDK = 35) DOUBLE PRECISION T(LDK), A(0:LDA), X, VAL, + BS, ONE, ZERO, LEFT, RIGHT, TOP REAL S INTRINSIC REAL EXTERNAL BS, FINDIT COMMON /KNOTS/ NKNOT, T, LEFT, RIGHT, TOP + /COEF/ A + /VALUE/ VAL + /INDEX/ INT C FIRST EXECUTABLE LINE ONE = 1.0D0 ZERO = 0.0D0 S = REAL (ZERO) CALL FINDIT (X, T, NKNOT, INT, LDK) IF (INT .EQ. 0) THEN IF (X .LT. LEFT) THEN S = REAL (ZERO) ELSE C A(0) IS CONSTANT TERM = S(T(1)) S = REAL(A(0) * (X-LEFT) / (T(1)-LEFT)) ENDIF ELSE IF (INT .EQ. NKNOT) THEN IF (X .GT. RIGHT) THEN S = REAL (ONE) ELSE S = REAL(TOP + (ONE - TOP) * (X-T(NKNOT)) / + (RIGHT - T(NKNOT))) ENDIF ELSE DO 100 I=0,NKNOT+1 S = S + REAL(A(I) * BS(X,I)) 100 CONTINUE ENDIF S = REAL(S - VAL) RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION BS(X,K) INTEGER NKNOT, INT, LDK, K PARAMETER ( LDK = 35) DOUBLE PRECISION T(LDK), Y, D, X, + ZERO, ONE, BS COMMON /KNOTS/ NKNOT, T + /INDEX/ INT DATA ONE /1.0D0/, + ZERO /0.0D0/ C Y IS A FUNCTION Y(D) = (3.0D0 - 2.0D0 * D) * (D**2) C FIRST EXECUTABLE LINE IF (K .EQ. 0) THEN BS = ONE ELSE IF (INT .LE. K-3) THEN BS = ZERO ELSE IF (INT .GE. K+1) THEN BS = ONE ELSE IF (K .EQ. 1) THEN D = (X - T(1)) / (T(2) - T(1)) BS = ONE - (ONE-D)**3 ELSE IF (K .EQ. 2) THEN D = (X - T(1)) / (T(3) - T(1)) BS = Y(D) ELSE IF (K .EQ. NKNOT) THEN D = (X - T(NKNOT-2)) / (T(NKNOT) - T(NKNOT-2)) BS = Y(D) ELSE IF (K .EQ. NKNOT+1) THEN D = (X - T(NKNOT-1)) / (T(NKNOT) - T(NKNOT-1)) BS = D**3 ELSE D = (X - T(K-2)) / (T(K+1) - T(K-2)) BS = Y(D) ENDIF RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION DS(X) INTEGER NKNOT, LDA, I, LDK, INT PARAMETER ( LDA = 37, + LDK = 35 ) DOUBLE PRECISION A(0:LDA),X, DBS, T(LDK), LEFT, RIGHT, + TOP, ZERO, THREE, ONE REAL DS EXTERNAL DBS, FINDIT INTRINSIC REAL COMMON /COEF/ A + /KNOTS/ NKNOT, T, LEFT, RIGHT,TOP + /INDEX/ INT ZERO = 0.0D0 ONE = 1.0D0 THREE = 3.0D0 DS = REAL (ZERO) CALL FINDIT (X, T, NKNOT, INT, LDK) IF (INT .EQ. 0) THEN IF (X .LT. LEFT) THEN DS = REAL(ZERO) ELSE DS = REAL(A(0) / (T(1)-LEFT)) ENDIF ELSE IF (INT .EQ. NKNOT) THEN IF (X .GT. RIGHT) THEN DS = REAL(ZERO) ELSE DS = REAL(A(NKNOT+1)*(ONE - TOP) / (RIGHT - T(NKNOT))) ENDIF ELSE DO 100 I=1,NKNOT+1 DS = DS + REAL (A(I) * DBS(X,I)) 100 CONTINUE ENDIF END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION DBS(X,K) INTEGER NKNOT, INT, LDK, K PARAMETER ( LDK = 35) DOUBLE PRECISION T(LDK), DY, D, X, + ZERO, ONE, DBS COMMON /KNOTS/ NKNOT, T + /INDEX/ INT DATA ONE /1.0D0/, + ZERO /0.0D0/ C Y IS A FUNCTION DY(D) = (ONE - D) * D * 6.0D0 C FIRST EXECUTABLE LINE IF (K .EQ. 0) THEN DBS = ZERO ELSE IF (INT .LE. K-3) THEN DBS = ZERO ELSE IF (INT .GE. K+1) THEN DBS = ZERO ELSE IF (K .EQ. 1) THEN D = (X - T(1)) / (T(2) - T(1)) DBS = 3.0D0*(ONE-D)**2 / (T(2)-T(1)) ELSE IF (K .EQ. 2) THEN D = (X - T(1)) / (T(3) - T(1)) DBS = DY(D) / (T(3) - T(1)) ELSE IF (K .EQ. NKNOT) THEN D = (X - T(NKNOT-2)) / (T(NKNOT) - T(NKNOT-2)) DBS = DY(D) / (T(NKNOT)-T(NKNOT-2)) ELSE IF (K .EQ. NKNOT+1) THEN D = (X - T(NKNOT-1)) / (T(NKNOT) - T(NKNOT-1)) DBS = 3.0D0*D**2 / (T(NKNOT)-T(NKNOT-1)) ELSE IF ((K .LE. NKNOT-2) .OR. (K .GE. 3)) THEN D = (X - T(K-2)) / (T(K+1) - T(K-2)) DBS = DY(D) / (T(K+1)-T(K-2)) ELSE WRITE (*,*) ' ERROR IN EVALUATING S; X NOT FOUND IN ANY', + ' INTERVAL.' ENDIF RETURN END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION DDS(X) INTEGER NKNOT, LDA, I, LDK, INT PARAMETER ( LDA = 37, + LDK = 35 ) DOUBLE PRECISION A(0:LDA), X, DDBS, T(LDK), LEFT, RIGHT, TOP, + ZERO REAL DDS EXTERNAL DDBS, FINDIT INTRINSIC REAL COMMON /COEF/ A + /KNOTS/ NKNOT, T, LEFT, RIGHT, TOP + /INDEX/ INT ZERO = 0.0D0 DDS = REAL (ZERO) CALL FINDIT (X, T, NKNOT, INT, LDK) IF ((INT .EQ. 0) .OR. (INT .EQ.NKNOT)) THEN DDS = REAL(ZERO) ELSE DO 100 I=1,NKNOT+1 DDS = DDS + REAL (A(I) * DDBS(X,I)) 100 CONTINUE ENDIF END ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION DDBS(X,K) INTEGER NKNOT, INT, LDK, K PARAMETER ( LDK = 35) DOUBLE PRECISION T(LDK), DDY, D, X, + ZERO, ONE, DDBS COMMON /KNOTS/ NKNOT, T + /INDEX/ INT DATA ONE /1.0D0/, + ZERO /0.0D0/ C Y IS A FUNCTION DDY(D) = (ONE - 2.0D0*D) * 6.0D0 C FIRST EXECUTABLE LINE IF (K .EQ. 0) THEN DDBS = ZERO ELSE IF (INT .LE. K-3) THEN DDBS = ZERO ELSE IF (INT .GE. K+1) THEN DDBS = ZERO ELSE IF (K .EQ. 1) THEN D = (X - T(1)) / (T(2) - T(1)) DDBS = -6.0D0*(ONE-D) / (T(2)-T(1))**2 ELSE IF (K .EQ. 2) THEN D = (X - T(1)) / (T(3) - T(1)) DDBS = DDY(D) / (T(3) - T(1))**2 ELSE IF (K .EQ. NKNOT) THEN D = (X - T(NKNOT-2)) / (T(NKNOT) - T(NKNOT-2)) DDBS = DDY(D) / (T(NKNOT)-T(NKNOT-2))**2 ELSE IF (K .EQ. NKNOT+1) THEN D = (X - T(NKNOT-1)) / (T(NKNOT) - T(NKNOT-1)) DDBS = 6.0D0*D / (T(NKNOT)-T(NKNOT-1))**2 ELSE D = (X - T(K-2)) / (T(K+1) - T(K-2)) DDBS = DDY(D) / (T(K+1)-T(K-2))**2 ENDIF RETURN END APPENDIX E: THE CALL STATEMENTS FROM SMOOTH 1: PROGRAM SMOOTH 91: CALL SDGEN 93: CALL GETDAT (X, NDATA) 94: CALL EPSSET (MAXFN, EPS, NSIG, RHO) 96: CALL CUMGEN (X, KNOTS, F, NDATA, NKNOT) 98: CALL SETWT (WEIGHT, NDATA) 106: CALL BSVALS (X, KNOTS, WEIGHT, NDATA, NKNOT, SVAL) 107: CALL CUMFVL (CUMF, WEIGHT, NDATA) 108: CALL VTPROF (SVAL, NDATA, NKNOT+2, LDX, SS) 109: CALL VMULFM (SVAL, CUMF, NDATA, NKNOT+2, 1, LDX, LDX, A, LDA, IER) 112: CALL LEQ1S (SS, NKNOT+2, A, 1, LDA, IJOB, IWK, WORK, IERLEQ) 113: CALL PLOT(X, KNOTS, NDATA, NKNOT, RUN) 115: CALL FIXSPL (FIX, RUN, RHO) 120: CALL ENDPTS 121: CALL FCOMP (KNOTS, F, NKNOT) 122: CALL GETSIZ (MAXSAM, SMPSZE, LDS, LDBT) 124: CALL GENSMD (F, NSIG, MAXFN, EPS, SAM, SMPSZE, WORK) 126: CALL GENBT (X, NDATA, SAM, SMPSZE) 129: CALL OUTPUT (X, F, NDATA, SMBT, BOOT, SMPSZE, 130: + MAXSAM, START, A) 134: SUBROUTINE SDGEN 153: SUBROUTINE GETDAT (X, N) 170: CALL GGEXN (DSEED, LAMBDA, N, X) 193: SUBROUTINE EPSSET (MAXFN, EPS, NSIG, RHO) 223: SUBROUTINE CUMGEN (X, X1, F, N, M) 259: CALL DETSTP (CRIT, START, STOPIT, N, LDK) 393: SUBROUTINE DETSTP (CRIT, START, STOP, N, LDK) 424: SUBROUTINE SETWT (WEIGHT, NDATA) 436: SUBROUTINE BSVALS (X, T, WEIGHT, NX, NK, SVAL) 447: CALL FINDIT (X(I), T, NK, INT, LDK) 455: SUBROUTINE CUMFVL (CUMF, WEIGHT, NDATA) 467: SUBROUTINE PLOT (X, T, NDATA, NABSI, RUN) 510: CALL UGETIO (3, 0, 8) 514: CALL USPLO(PTWORK(0,1),PTWORK(0,2),1+NINT,NINT+1,1,1,ITITLE, 515: + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) 518: CALL USPLO(PTWORK(0,1),PTWORK(0,3),1+NINT,NINT+1,1,1,ITITLE, 519: + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) 522: CALL USPLO(PTWORK(0,1),PTWORK(0,4),1+NINT,NINT+1,1,1,ITITLE, 523: + NTITLE,IXLABL,NXLABL,IYLABL,NYLABL,RANGE,ICHAR,IOPT,IER) 529: SUBROUTINE FIXSPL (FIX, RUN, LAMBDA) 595: CALL EXPSM (ADDX, LAMBDA) 601: CALL EXPSM (X(NX), LAMBDA) 607: CALL EXPSM (X(1), LAMBDA) 619: SUBROUTINE EXPSM (PIVOT, LAMBDA) 642: SUBROUTINE ENDPTS 674: SUBROUTINE FCOMP (X, F, NF) 688: SUBROUTINE GETSIZ (NBOOTS, SZBOOT, LDS, LDBT) 702: SUBROUTINE GENSMD (F, NSIG, ITMAX, EPS, SAMP, K, WORK) 725: CALL GGUBS (DSEED, K, WORK) 729: CALL FINDIT (SAMP(I), F, NABSI, INTERV, LDK) 746: CALL ZFALSE (S, EPS, NSIG, T1, T2, SAMP(I), ITMAX, IER) 756: SUBROUTINE GENBT (X, N, SAMP, M) 771: CALL GGUD (DSEED, N, M, ISAMP) 793: SUBROUTINE OUTPUT(X, F, NX, SMTHBT, REGBT, SZBT, 794: + NBT, START, A) 871: SUBROUTINE FINDIT (VAL, Y, M, IT, LDY) 922: CALL FINDIT (X, T, NKNOT, INT, LDK) 999: CALL FINDIT (X, T, NKNOT, INT, LDK) 1073: CALL FINDIT (X, T, NKNOT, INT, LDK) The following subroutines were called from IMSL: VTPROF VMULFM LEQ1S GGEXN UGETIO USPLO GGUBS ZFALSE GGUD