PROGRAM PICTURE c--- This collection of routines provides a picture of the data c set by displaying: c c (a) original data (sorted) c (b) a stem-and-leaf plot c (c) a table of the quantiles, spreads, and ratios c (d) a plot of the midspreads vs quantiles c (e) a hanging histogram of the data c***************************************************************************** integer maxin,err,iunit,iw(1000) integer p(130),pmax,pmin,outptr,maxptr,ounit integer linset,chrset logical xtrems real epsi,maxint,mhat,shat real xdata(1000),data(1000),sorty(1000) real wx(10),wy(10),x(1000),y(1000),yhat(1000) real drr(1000),gsp(9),ylv(15,2),clv(8),cclv(1000,8) real d(15) character*32 filename,outputfile data gsp /1.,1.349,2.301,3.068,3.726,4.308,4.836,5.320,5.320/ common/chrbuf/ p,pmax,pmin,outptr,maxptr,ounit common/numbrs/epsi,maxint c--- initialize the EDA routines iunit = 5 ounit = 20 pmin = 2 pmax = 80 epsi = 1.e-7 maxin = 32000 err = 0 CALL CINIT (OUNIT,PMIN,PMAX,EPSI,MAXIN,ERR) c--- read in the data print *, 'Input data file: ' read(iunit,'(a)') filename print *, 'Output data file: ' read (iunit,'(a)') outputfile open (unit = iunit,status='old',file=filename) open (unit=ounit,status='new',file=outputfile) n = 0 do 5 i=1,1000 read (5,*,end=10,err=10) data(i) n = n + 1 5 continue c--- sort the data and obtain the letter values 10 CALL LVALS (DATA,N,D,YLV,NLV,XDATA,ERR) c--- print out the sorted data write (ounit,8) filename 8 format (t25,'FILE = ',a) write (ounit,7) write (ounit,9) n 9 format(t25,'SORTED DATA, N = ',i3) write (ounit,7) write (ounit,11) (xdata(i),i=1,n) 11 format (5(1x,f10.3)) c--- call the stem and leaf routine xtrems = .false. write (ounit,7) CALL STMLF (XDATA,N,SORTY,IW,XTREMS,ERR) c--- obtain corrections to letter values for PSEUDOSPREADS if(nlv.gt.8) nlv = 8 ! throw away other letter values c--- if n > 100 then use the approx. corrections to lv if (n.gt.100) then xn=float(n) cclv(n,1)=1.-(2./xn)*(1.-.86*((xn/2.)-(n/2))) cclv(n,2)=1.71-(3./xn)*(1.-.8*((xn/4.)-(n/4))) cclv(n,3)=2.27-(5./xn)*(1.-.77*((xn/8.)-(n/8))) cclv(n,4)=2.76-(10./xn)*(1.-.76*((xn/16.)-(n/16))) cclv(n,5)=3.19-(15./xn)*(1.-.76*((xn/32.)-(n/32))) cclv(n,6)=3.58-(30./xn)*(1.-.75*((xn/64.)-(n/64))) cclv(n,7)=3.94-(55./xn)*(1.-.75*((xn/128.)-(n/128))) do i=1,nlv-1 cclv(n,i)=1.349*cclv(n,i) enddo else open(unit=8,file='corr.dat',status='old') do i=1,100 read (8,*) id,(cclv(i,j),j=1,7) cclv(i,8) = 1.0 enddo rewind (unit=8) endif c--- obtain corrections to spreads as a function of n do 40 i=2,nlv clv(i) = cclv(n,i-1) 40 continue c--- print out the leter values, sreads, and ratios write (ounit,71) write (ounit,18) 18 format(t20,'LETTER VALUE MATRIX') write (ounit,71) 6 format('1') 7 format(//) 71 format(/) CALL LVPRNT (NLV,D,YLV,CLV,ERR) c--- print out a plot of the midsummary vs z**2 linset = 15 chrset = 10 xmax= -1.e8 ymax = -1.e8 do i=1,nlv x(i) = (gsp(i)/2.)**2 y(i) = (ylv(i,1)+ylv(i,2))/2. ymax = amax1(ymax,y(i)) xmax = amax1(xmax,x(i)) enddo write (ounit,6) write (ounit,7) write (ounit,8) filename write (ounit,71) write (ounit,15) 15 format(t20,'PLOT OF MIDSUMMARY VS ZEE**2') write (ounit,71) ymin=ymax xmin=xmax CALL PLOT (Y,X,NLV,WY,WX,LINSET,CHRSET, + XMIN,XMAX,YMIN,YMAX,ERR) c--- print out a plot of the gaussian spread vs z**2 linset = 15 chrset = 10 xmax= -1.e8 ymax = -1.e8 do i=2,nlv x(i-1) = (gsp(i)/2.)**2 y(i-1) = abs(ylv(i,2)-ylv(i,1))/clv(i) ymax = amax1(ymax,y(i-1)) xmax = amax1(xmax,x(i-1)) enddo write (ounit,7) write (ounit,14) 14 format(t20,'PLOT OF PSEUDOSIGMA VS ZEE**2') write (ounit,71) ymin=ymax xmin=xmax CALL PLOT (Y,X,NLV-1,WY,WX,LINSET,CHRSET, + XMIN,XMAX,YMIN,YMAX,ERR) c--- choose reasonable bins for data fsig = abs(ylv(2,2)-ylv(2,1)) bw = 2*fsig/(n**0.333) nbins = abs(ylv(nlv,2)-ylv(nlv,1))/bw c--- bin the data xlow = ylv(nlv,1) do i=1,nbins+1 y(i)=0. x(i) = xlow + (i-1)*bw enddo do i=1,n do j=1,nbins if (xdata(i).ge.x(j).and.xdata(i).lt.x(j+1)) + y(j) = y(j)+1 enddo enddo c--- call the rootogram routine write (ounit,6) write (ounit,8) filename write (ounit,7) write (ounit,16) 16 format(t20,'ROOTOGRAM OF FITTED GAUSSIAN RESIDUALS') write (ounit,7) CALL RGCOMP (X,Y,NBINS,MU,SIGMA,YHAT,DRR,MHAT,SHAT,ERR) CALL RGPRNT (Y,NBINS,YHAT,DRR,ERR) c--- print out the fitted parameters write (ounit,7) write (ounit,12) mhat write (ounit,13) shat 12 format(' Fitted mean = ',f10.3) 13 format(' Fitted standard deviation = ',f10.3) stop end BLOCK DATA c--a group of utility routines which are required for the c other EDA programs c--- chars contains the symbols of the standard fortran character c set, and cha - chpt are the corresponding indices into chars. c putchr is the primary user of this translation vector. common /chario/ chars,cmax, + cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk, + chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv, + chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6, + ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh, + chlpar,chrpar,chcoma,chpt integer chars(46),cmax integer cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk integer chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv integer chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6 integer ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh integer chlpar,chrpar,chcoma,chpt data chars( 1),chars( 2),chars( 3),chars( 4) /1HA,1HB,1HC,1HD/ data chars( 5),chars( 6),chars( 7),chars( 8) /1HE,1HF,1HG,1HH/ data chars( 9),chars(10),chars(11),chars(12) /1HI,1HJ,1HK,1HL/ data chars(13),chars(14),chars(15),chars(16) /1HM,1HN,1HO,1HP/ data chars(17),chars(18),chars(19),chars(20) /1HQ,1HR,1HS,1HT/ data chars(21),chars(22),chars(23),chars(24) /1HU,1HV,1HW,1HX/ data chars(25),chars(26),chars(27),chars(28) /1HY,1HZ,1H0,1H1/ data chars(29),chars(30),chars(31),chars(32) /1H2,1H3,1H4,1H5/ data chars(33),chars(34),chars(35),chars(36) /1H6,1H7,1H8,1H9/ data chars(37),chars(38),chars(39),chars(40) /1H ,1H=,1H+,1H-/ data chars(41),chars(42),chars(43),chars(44) /1H*,1H/,1H(,1H)/ data chars(45),chars(46) /1H,,1H./ data cmax/46/ data cha,chb,chc,chd,che,chf /1,2,3,4,5,6/ data chg,chh,chi,chj,chk,chl /7,8,9,10,11,12/ data chm,chn,cho,chp,chq,chr /13,14,15,16,17,18/ data chs,cht,chu,chv,chw,chx /19,20,21,22,23,24/ data chy,chz,ch0,ch1,ch2,ch3 /25,26,27,28,29,30/ data ch4,ch5,ch6,ch7,ch8,ch9 /31,32,33,34,35,36/ data chbl,cheq,chplus,chmin /37,38,39,40/ data chstar,chslsh,chlpar,chrpar /41,42,43,44/ data chcoma,chpt /45,46/ end subroutine cinit (iounit,ipmin,ipmax,iepsi,imaxin,err) integer iounit,ipmin,ipmax,imaxin,err real iepsi c--- initialization, to be called at start of any main program c which calls one of the EDA subroutines c--- IOUNIT is the number of the unit wo thich output is directed c IPMIN is the left margin c IPMAX is the right margin c IEPSI is the machine-realted epsilon c IMAXIN is the maximum permitted integer value c ERR is the (usual) error flag, to indicate whether c the routine executed successfully integer p(130), pmax,pmin,outptr,maxptr,ounit real epsi,maxint common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit common /numbrs/ epsi,maxint integer blank,i data blank /1H / err = 6 if(ipmin .lt. 1) goto 999 if(ipmax .gt. 130) goto 999 if(ipmax .le. ipmin) goto 999 err = 7 if((1.0 + iepsi) .le. 1.0) goto 999 err = 0 ounit = iounit pmin = ipmin outptr = ipmin maxptr = ipmin pmax = ipmax epsi = iepsi maxint = float(imaxin) do 50 i = 1,130 p(i) = blank 50 continue 999 return end subroutine putchr (posn,char,err) integer posn,char,err c--- Place the character CHAR at position POSN in the c output line P. If POSN = 0, place CHAR in the c next available position in P. MAXPTR is to be initialized c to PMIN, and PRINT must reset it. common /chario/ chars,cmax, + cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk, + chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv, + chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6, + ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh, + chlpar,chrpar,chcoma,chpt integer chars(46),cmax integer cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk integer chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv integer chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6 integer ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh integer chlpar,chrpar,chcoma,chpt integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit if(char .gt. 0 .and. char .le. cmax) goto 10 err = 4 return 10 if(posn .ne. 0) outptr = max0(pmin,posn) outptr = min0(outptr,pmax) p(outptr) = chars(char) maxptr = max0(maxptr,outptr) outptr = outptr + 1 return end integer function wdthof(i) integer i c--- find the number of characters needed to print I integer ia,iq,nd ia = iabs(i) nd = 1 if(i.lt. 0) nd = 2 10 iq = ia/10 if(iq.eq.0) goto 20 ia = iq nd = nd + 1 goto 10 20 wdthof = nd return end subroutine putnum (posn,n,w,err) integer posn,n,w,err c--- Place the character representation of the integer n c right justified in a field w spaces wide starting at position c posn in the output line p. c The variables ip,inum, and iw are internal versions of c posn, n, and w. We proceed by extracting the digits c of n starting with the low-order digit, and stacking them in c dstk. (nd counts the digits). Once we have collected all the digits c (and know that w spaces are sufficient), we skip over any unneeded c spaces, put out a minus sugn if needed, and then put out c the digits, starting with the high-order one. c This routine calls putchr and depends on having digits c 0 through 9 in consecutive elements of chars in the c common block chario, starting at ch0 = 27. It also c assumes that the minus sign is at chmin = 40 in chars. integer chd,ch0,chmin,dstk(20),inum,ip,iq,iw,nd integer p(130), pmax,pmin,outptr,maxptr,ounit common/chrbuf/ p,pmax,pmin,outptr,maxptr,ounit data ch0,chmin/27,40/ iw = w if(n .lt. 0) iw = iw - 1 inum = iabs(n) c--- extract and stack the digits of inum, checking to see that c n fits in w spaces nd = 1 10 iq = inum/10 dstk(nd) = inum - iq*10 if(nd .le. 20 .and. nd .le. iw) goto 20 err = 2 goto 999 20 if(iq .eq. 0) goto 30 inum = iq nd = nd + 1 goto 10 c--- unstack the digits from dstk and put them out. c Note that when n is negative, a minus sign must be c inserted in the space before the first digit. Decreasing c iw by 1 in the initialization has provided a space for the c minus sign 30 ip = posn if(ip .eq. 0) ip = outptr ip = ip + iw - nd if(n .ge. 0) goto 40 call putchr (ip,chmin,err) ip = ip + 1 40 chd = ch0 + dstk(nd) call putchr (ip,chd,err) if(nd .eq. 1) goto 50 nd = nd - 1 ip = ip + 1 goto 40 50 continue 999 return end subroutine print c--- Print the output line p on unit ounit (maxptr c indicates the rightmost position which has been used c in this line). Then reset p to spaces and maxptr and outptr to c pmin. integer p(130),pmax,pmin,outptr,macptr,ounit common /chrbuf/p,pmax,pmin,outptr,maxptr,ounit integer blank,i data blank/1H / write (ounit,10) (p(i),i= 1, maxptr) 10 format(1x,130a1) do 20 i=1,maxptr p(i) = blank 20 continue outptr = pmin maxptr = pmin return end SUBROUTINE QSORT8(N,X) * Sorting program that uses a quicksort algorithm REAL*8 X(N) REAL*8 KEY, KL, KR, KM, TEMP INTEGER*4 L, R, M, LSTACK(50), RSTACK(50), SP INTEGER*4 NSTOP LOGICAL*1 MGTL, LGTR, RGTM nstop = 15 IF(N.LE.NSTOP) GOTO 100 SP = 0 SP = SP + 1 LSTACK(SP) = 1 RSTACK(SP) = N * Sort a subrecord off the stack * Set KEY = median of X(L), X(M), X(R) 1 L = LSTACK(SP) R = RSTACK(SP) SP = SP - 1 M = (L + R) / 2 KL = X(L) KM = X(M) KR = X(R) MGTL = KM .GT. KL RGTM = KR .GT. KM LGTR = KL .GT. KR IF(MGTL .EQV. RGTM) THEN IF(MGTL .EQV. LGTR) THEN KEY = KR ELSE KEY = KL ENDIF ELSE KEY = KM ENDIF I = L J = R * Find a big record on the left 10 IF(X(I).GE.KEY) GOTO 11 I = I + 1 GOTO 10 11 CONTINUE * Find a small record on the right 20 IF(X(J).LE.KEY) GOTO 21 J = J - 1 GOTO 20 21 CONTINUE IF(I.GE.J) GOTO 2 * Exchange records TEMP = X(I) X(I) = X(J) X(J) = TEMP I = I + 1 J = J - 1 GOTO 10 * Subfile is partitioned into two halves, left .le. right * Push the two halves on the stack 2 IF(J-L+1 .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = L RSTACK(SP) = J ENDIF IF(R-J .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = J+1 RSTACK(SP) = R ENDIF * Anything left to process? IF(SP.GT.0) GOTO 1 * Sorting routine that sorts the N elements of double precision * array X by straight insertion between previously sorted numbers 100 DO 110 J = N-1,1,-1 K = J DO 120 I = J+1,N IF(X(J).LE.X(I)) GOTO 121 120 K = I 121 CONTINUE IF(K.EQ.J) GOTO 110 TEMP = X(J) DO 130 I = J+1,K 130 X(I-1) = X(I) X(K) = TEMP 110 CONTINUE RETURN END SUBROUTINE QSORT4(N,X) * Sorting program that uses a quicksort algorithm REAL*4 X(N) REAL*4 KEY, KL, KR, KM, TEMP INTEGER*4 L, R, M, LSTACK(50), RSTACK(50), SP INTEGER*4 NSTOP LOGICAL*1 MGTL, LGTR, RGTM nstop = 15 IF(N.LE.NSTOP) GOTO 100 SP = 0 SP = SP + 1 LSTACK(SP) = 1 RSTACK(SP) = N * Sort a subrecord off the stack * Set KEY = median of X(L), X(M), X(R) 1 L = LSTACK(SP) R = RSTACK(SP) SP = SP - 1 M = (L + R) / 2 KL = X(L) KM = X(M) KR = X(R) MGTL = KM .GT. KL RGTM = KR .GT. KM LGTR = KL .GT. KR IF(MGTL .EQV. RGTM) THEN IF(MGTL .EQV. LGTR) THEN KEY = KR ELSE KEY = KL ENDIF ELSE KEY = KM ENDIF I = L J = R * Find a big record on the left 10 IF(X(I).GE.KEY) GOTO 11 I = I + 1 GOTO 10 11 CONTINUE * Find a small record on the right 20 IF(X(J).LE.KEY) GOTO 21 J = J - 1 GOTO 20 21 CONTINUE IF(I.GE.J) GOTO 2 * Exchange records TEMP = X(I) X(I) = X(J) X(J) = TEMP I = I + 1 J = J - 1 GOTO 10 * Subfile is partitioned into two halves, left .le. right * Push the two halves on the stack 2 IF(J-L+1 .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = L RSTACK(SP) = J ENDIF IF(R-J .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = J+1 RSTACK(SP) = R ENDIF * Anything left to process? IF(SP.GT.0) GOTO 1 * Sorting routine that sorts the N elements of single precision * array X by straight insertion between previously sorted numbers 100 DO 110 J = N-1,1,-1 K = J DO 120 I = J+1,N IF(X(J).LE.X(I)) GOTO 121 120 K = I 121 CONTINUE IF(K.EQ.J) GOTO 110 TEMP = X(J) DO 130 I = J+1,K 130 X(I-1) = X(I) X(K) = TEMP 110 CONTINUE RETURN END SUBROUTINE QSORT2(N,X,Y) * Sorting program that uses a quicksort algorithm, sorts both x and y by x REAL*4 X(N), Y(N) REAL*4 KEY, KL, KR, KM, TEMP INTEGER*4 L, R, M, LSTACK(50), RSTACK(50), SP INTEGER*4 NSTOP LOGICAL*1 MGTL, LGTR, RGTM nstop = 15 IF(N.LE.NSTOP) GOTO 100 SP = 0 SP = SP + 1 LSTACK(SP) = 1 RSTACK(SP) = N * Sort a subrecord off the stack * Set KEY = median of X(L), X(M), X(R) 1 L = LSTACK(SP) R = RSTACK(SP) SP = SP - 1 M = (L + R) / 2 KL = X(L) KM = X(M) KR = X(R) MGTL = KM .GT. KL RGTM = KR .GT. KM LGTR = KL .GT. KR IF(MGTL .EQV. RGTM) THEN IF(MGTL .EQV. LGTR) THEN KEY = KR ELSE KEY = KL ENDIF ELSE KEY = KM ENDIF I = L J = R * Find a big record on the left 10 IF(X(I).GE.KEY) GOTO 11 I = I + 1 GOTO 10 11 CONTINUE * Find a small record on the right 20 IF(X(J).LE.KEY) GOTO 21 J = J - 1 GOTO 20 21 CONTINUE IF(I.GE.J) GOTO 2 * Exchange records TEMP = X(I) X(I) = X(J) X(J) = TEMP TEMP = Y(I) Y(I) = Y(J) Y(J) = TEMP I = I + 1 J = J - 1 GOTO 10 * Subfile is partitioned into two halves, left .le. right * Push the two halves on the stack 2 IF(J-L+1 .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = L RSTACK(SP) = J ENDIF IF(R-J .GT. NSTOP) THEN SP = SP + 1 LSTACK(SP) = J+1 RSTACK(SP) = R ENDIF * Anything left to process? IF(SP.GT.0) GOTO 1 * Sorting routine that sorts the N elements of single precision * array X by straight insertion between previously sorted numbers 100 DO 110 J = N-1,1,-1 K = J DO 120 I = J+1,N IF(X(J).LE.X(I)) GOTO 121 120 K = I 121 CONTINUE IF(K.EQ.J) GOTO 110 TEMPX = X(J) TEMPY = Y(J) DO 130 I = J+1,K X(I-1) = X(I) 130 Y(I-1) = Y(I) X(K) = TEMPX Y(K) = TEMPY 110 CONTINUE RETURN END subroutine yinfo (y,n,med,hl,hh,adjl,adjh,iadjl,iadjh,step,err) c--- Gets general information about Y(). Useful for plot scaling. c Sorts Y() and returns it sorted. Also return: c med = median c hl = low hinge c hh = hi hinge c adjl = low adjacent value c adjh = hi adjacent value c iadjl = location of adjl c iadjh = location of adjh integer n,iadjl,iadjh,err real y(n),med,hl,hh,adjl,adjh,step real hfence,lfence integer j,k,temp1,temp2 call qsort4(n,y) k = n j = (k/2)+1 temp1 = n+1-j med = (y(j) + y(temp1))/2.0 k = (k+1)/2 j = (k/2) + 1 temp1 = k+1 - j hl = (y(j) + y(temp1))/2.0 temp1 = n-k+j temp2 = n+1-j hh = (y(temp1)+y(temp2))/2.0 step = (hh - hl)*1.5 hfence = hh + step lfence = hl - step c--- find adjacent values iadjl = 0 20 iadjl = iadjl + 1 if (y(iadjl) .le. lfence) goto 20 adjl = y(iadjl) iadjh = n+1 30 iadjh = iadjh - 1 if(y(iadjh) .ge. hfence) goto 30 adjh = y(iadjh) 999 return end subroutine nposw (hi,lo,nionos,nn,maxp,mzero,ptotl,fract, + unit,npw,err) c--- Find a nice (i.e. simple) data units value to assign to one plot c position in one dimension of a plot. A plot position is typically c one character position horizontally, or one line vertically. c On entry: c hi,lo the high and low edges of the data range to be plotted c nionos a vector of length nn containing nice mantissas for the c plot unit c maxp maximum number of plot positions allowed in this dimension c of the plot c mzero is .true. if a position labelled -0 is allowed in this c dimension, .false. otherwise c c On exit: c ptotl holds the total number of plot positions to be used c in this dimension. c fract is the mantissa of the nice position width. It is c selected from nionos c unit is an integer power of 10 such that npw = fract*unit c npw is the nice position width. One plot position width c will represent a data-space distance of npw. integer nn,maxp,ptotl,err real hi,lo,nionos(nn),fract,unit,npw logical mzero integer floor,intfn integer i real aprxw if (maxp .gt. 0) goto 5 err = 8 goto 999 5 aprxw = (hi - lo)/float(maxp) if(aprxw .gt. 0.0) goto 10 err = 9 goto 999 10 unit = 10.0 **floor(alog10(aprxw)) fract = aprxw/unit do 20 i=1,nn if(fract .le. nionos(i)) goto 30 20 continue 30 fract = nionos(i) npw = fract * unit ptotl = intfn(hi/npw,err) - intfn(lo/npw,err) + 1 if(err .ne. 0) goto 999 c--- if minus zero position possible and sgn(hi) .ne. sgn(lo) allow it if(mzero .and. (hi*lo.lt.0.0.or.hi.eq.0.0)) ptotl = ptotl+1 c--- ptotl positions required with this width -- few enough? if(ptotl .le. maxp) goto 999 c--- too many positions needed, so bump npw up one nice number i = i+1 if(i.le. nn) goto 30 i = 1 unit = unit * 10.0 goto 30 999 return end integer function intfn (x,err) c--- find the integer equal to or next closer to zero than x c--- check to see that x is not too large to fit in an integer variable real x integer err common /numbrs/ epsi,maxint real epsi,maxint if(abs(x) .le. maxint) goto 10 c--- x is too large in magnitude to fit into an integer, so return c the largest legal integer and set the error flag err = 3 intfn = ifix(sign(maxint,x)) goto 999 10 intfn = int((1.0 + epsi) * x) 999 return end integer function floor (y) real y c--- find floor(y) the largest integer not exceeding y floor = int(y) if(y .lt. 0.0 .and. y .ne. float(floor)) floor = floor - 1 return end real function median (y,n) c-- find the median of the sorted values y() integer n real y(n) integer mptr,mpt2 mptr = (n/2) + 1 mpt2 = n-mptr+1 median = (y(mptr) + y(mpt2))/2.0 return end real function gau(z) c--- This function calculates the value of the standard c Gaussian cumulative distribution function at z. c The algorith uses approximations given in c Mathematics of Computation, v. 31 (1977) p. 214 real p,pi,x,z x = abs(z) if(x.gt.5.5) goto 10 p = exp(-((83.0 * x + 351.0) * x + 562.0) * x / + (703.0 + 165.0 * x)) goto 20 10 pi = 4.0 * atan(1.0) p = sqrt(2.0/pi)*exp(-(x*x/2.0 + 0.94/(x*x)))/X c--- the approximations yield values of the half-normal tails c so we translate that into the value of the Gaussian c.d.f. c and allow for the sign of z 20 gau = p/2.0 if(z.gt.0.0) gau = 1.0 - gau return end subroutine lvals (y,n,d,ylv,nlv,sorty,err) integer n,nlv,err real y(n),d(15),ylv(15,2),sorty(n) c--- For the batch of values in y, find the selected quantiles known c as the letter values. Upon exit, ylv contains the letter values, c d contains the corresponding depths, and nlv is the number of pairs c of letter values. Specifically, ylv(1,1) and ylv(1,2) are both c set equal to the median, whose depth d(1), is (n+1)/2. The rest of c the letter values come in pairs and are stored in ylv in order from c the hinges out to the extremes. Thus ylv(2,1) and ylv(2,2) are the c lower hinge and upper hinge, respectively, and ylv(nlv,1) and ylv(nlv,2) c are the lower extreme (minimum) and upper extremem (maximum), respectively. c--- local variables integer i,j,k,pt1,pt2 if((n .gt. 3) .and. (n .le. 24576)) goto 10 nlv = 0 err = 21 goto 999 c--- sort y into sorty 10 continue do 15 i=1,n sorty(i) = y(i) 15 continue call qsort4(n,sorty) c--- handle median separately because it is not a pair of lv d(1) = float(n+1)/2.0 j = (n/2) + 1 pt2 = n + 1 - j ylv(1,1) = (sorty(j) + sorty(pt2))/2.0 ylv(1,2) = ylv(1,1) k=n i=2 20 k= (k+1)/2 j = (k/2) + 1 d(i) = float(k+1)/2.0 pt2 = k + 1 - j ylv(i,1) = (sorty(j) + sorty(pt2))/2.0 pt1 = n-k+j pt2 = n+1-j ylv(i,2) = (sorty(pt1) + sorty(pt2))/2.0 i = i + 1 if(d(i-1) .gt. 2.0) goto 20 nlv = i d(i) = 1.0 ylv(i,1) = sorty(1) ylv(i,2) = sorty(n) 999 return end subroutine lvprnt (nlv,d,ylv,clv,err) integer nlv,err real d(15),ylv(15,2),clv(8) c--- prints a letter value display. c The nlv pairs of letter values are in ylv -- ylv(1,1) is the c lower letter value in the pair and ylv(i,2) is the upper c letter value, with the exception that ylv(1,1) and ylv (1,2) c are both equal to the median. The vector d contains the c corresponding depths. integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/p,pmax,pmin,outptr,maxptr,ounit c--- local variables integer i,n,tags(14) real mid,spr data tags(1),tags(2),tags(3),tags(4) /1HM,1HH,1HE,1HD/ data tags(5),tags(6),tags(7),tags(8) /1HC,1HB,1HA,1HZ/ data tags(9),tags(10),tags(11),tags(12) /1HY,1HX,1HW,1HV/ data tags(13),tags(14) /1HU,1HT/ if((nlv .ge. 3) .and. (nlv .le. 15)) goto 10 err = 22 goto 999 10 if(pmax .ge. 64) goto 20 err = 23 goto 999 20 write (ounit,1001) 1001 format(5x,'DEPTH',7X,'LOWER',8X,'UPPER',11X,'MID', + 8X,'SPREAD',8x,'PSEUDO') c--- recover n from D(1), the depth of the median n = int(2.0*d(1)) - 1 write (ounit,1002) n 1002 format(1x,'N=',i5) c--- write line containing median (and first mid) write (ounit,1003) d(1),ylv(1,1),ylv(1,1) 1003 format(1x,'M ',F7.1,8X,F10.3,13X,F10.3) n = nlv - 1 do 30 i = 2,n mid = (ylv(i,1) + ylv(i,2))/2.0 spr = ylv(i,2) - ylv(i,1) pseu = spr/clv(i) write (ounit,1004) tags(i),d(i),ylv(i,1), + ylv(i,2),mid,spr,pseu 1004 format(1x,a1,1x,f7.1,3x,f10.3,3x,f10.3,5x, + f10.3,3x,f10.3,3x,f10.3) 30 continue mid = (ylv(nlv,1) + ylv(nlv,2))/2.0 spr = ylv(nlv,2) - ylv(nlv,1) pseu = spr/clv(nlv) write (ounit,1005) ylv(nlv,1),ylv(nlv,2),mid,spr,pseu 1005 format (7x,'1',5x,f10.3,3x,f10.3,5x, + f10.3,3x,f10.3,3x,f10.3/) 999 return end subroutine stmlf (y,n,sorty,iw,xtrems,err) integer n,iw(n),err real y(n),sorty(n) logical xtrems c--- produce a stem-and-leaf display of the data in y(n) c--- IW is an iteger work array. Sorty() is a real work array. c xtrems is a logical flag, .true. if scaling to extremes. c Otherwise scales to fences. c--- common blocks and variables for output integer p(130),pmax,pmin,outptr,maxptr,ounit real epsi,maxint common/chrbuf/p,pmax,pmin,outptr,maxptr,ounit common/numbrs/epsi,maxint c--- functions integer intfn,floor c--- calls subroutines depthp, nposw, outlyp, putchr, putnum, c sltitl, stemp, yinfo real med,hl,hh,adjl,adjh,step,unit,nionos(4),npw integer i,slbrk,pltwid,rank,iadjl,iadjh,nlins integer nlmax,linwid integer low,hi,cut,stem,pt1,pt2,j,spacnt,leaf,nn,chstar logical negnow,medyet c--- data definitions data chstar/41/ data nionos(1),nionos(2),nionos(3),nionos(4)/1.0,2.0,5.0,10.0/ data nn/4/ if(n .ge. 2) goto 5 err = 11 goto 999 c--- setup: find width of plotting region, stem-leaf break position, etc. 5 slbrk = pmin +11 pltwid = pmax - slbrk - 2 if(pltwid .gt. 5) goto 10 err = 13 goto 999 c--- find the best scale for the plot c--- sort y in sorty and get summary information 10 do 20 i=1,n sorty(i) = y(i) 20 continue call yinfo (sorty,n,med,hl,hh,adjl,adjh,iadjl,iadjh,step,err) if(err .ne. 0) goto 999 c--- find nice line width for plot c--- if adjacent values equal or user demands it, fake the adjacent c values to be the extremes if((adjh .gt. adjl) .and. .not. xtrems) goto 25 iadjl = 1 iadjh = n adjl = sorty(iadjl) adjh = sorty(iadjh) c-- note that these used to be y() not sorty() 25 nlmax = intfn (10.0*alog10(float(iadjh - iadjl + 1)), err) if(adjh .gt. adjl ) goto 27 c--- even if all values are equal we can produce a display adjh = adjl + 1.0 nlmax = 1 27 call nposw (adjh,adjl,nionos,nn,nlmax,.true.,nlins,fract, + unit,npw,err) if(err .ne. 0) goto 999 c--- rescale everything according to unit. Hereafter everything is integer c and data are of the form ss...sl(.) c note that intfn performs epsil?n adjustments for correct rounding, c and checks that the real number is not too large for an integer varibale. do 30 i=1,n iw(i) = intfn(sorty(i)/unit,err) 30 continue if(err. ne. 0) goto 999 if(fract .eq. 10.0) goto 40 c--- if all leaves are zero, we should be in one-line-per-stem format do 35 i=iadjl,iadjh if(mod(iw(i),10) .ne. 0) goto 40 35 continue fract = 10.0 npw = fract * unit nlins = intfn (adjh/npw,err) - intfn (adjl/npw,err) + 1 if(adjh*adjl .lt. 0.0 .or. adjh .eq. 0.0) nlins = nlins + 1 40 low = iw(iadjl) hi = iw(iadjh) c--- linewidth now is nicewidth/unit = fract linwid = intfn (fract,err) call sltitl (unit,err) if(err .ne. 0) goto 999 c--- print values below low adjacent value on "lo" stem rank =iadjl - 1 if(iadjl .eq. 1) goto 50 call outlyp (iw,n,1,rank,.false.,slbrk,err) if(err .ne.0) goto 999 c--- initialize for main part of display c initial setting are to line before first one printed 50 cut = floor ((1.0 + epsi)*float(low)/float(linwid))*linwid negnow = .true. stem = cut if(low .lt. 0) goto 60 c--- first stem positive negnow = .false. stem = cut - linwid 60 medyet = .false. c--- Two pointers are used. PT1 counts first for depths, PT2 follows c for leaf printing. Both are initialized one point early pt1 = iadjl pt2 = pt1 c--- main loop. for each line do 120 j=1,nlins c--- variable uses: c CUT : first number on next line of positive stems c but last number on current line of negative stems c STEM: inner (near zero) edge of current line c SPACNT: counts spaces used on this line c--- step to next line cut = cut + linwid if(stem .ne. 0 .or. .not. negnow) goto 70 negnow = .false. goto 80 70 stem = stem + linwid c--- newline: initialie count of spaces used 80 spacnt = 0 c--- find and print depth call depthp (sorty,iw,n,pt1,pt2,cut,iadjh,hi,rank, + medyet,slbrk,err) if(err .ne. 0) goto 999 c--- print stem label call stemp (stem,linwid,negnow,slbrk,err) if(err .ne. 0) goto 999 c--- find and print leaves if(pt1 .eq. pt2) goto 110 90 leaf = iabs(iw(pt2) - (stem/10)*10) call putnum(0,leaf,1,err) spacnt = spacnt + 1 if(spacnt .lt. pltwid) goto 100 c--- line overflows past right edge, mark with * call putchr (0,chstar,err) if(err .ne. 0) goto 999 pt2 = pt1 goto 110 100 pt2 = pt2 + 1 if(pt2 .lt. pt1) goto 90 c--- end line 110 call print c--- continue loop until we run out of numbers to plot 120 continue c--- print values above hi adjacent value on "hi" stem if(pt1 .gt. n) goto 990 call outlyp (iw,n,pt1,n,.true.,slbrk,err) 990 write (ounit,5990) 5990 format(1x) 999 return end subroutine outlyp (iw,n,from,to,hiend,slbrk,err) logical hiend integer n,iw(n),from,to,slbrk,err c--- print the lo or hi stem for a stem-and-leaf display. c The logical variable HIEND is .true. if we are to print c the high stem, .false. if the lo stem is to be printed. c IW() contains n sorted and scaled data values. Each has the c form ss...sl, where th ones digit is the leaf. c FROM, TO are pointers into IW() delimiting the values to be c placed on the HI or LO stem. SLBRK is the character position c on the page of the blank column between stems and leaves. integer p(130),pmax,pmin,outptr,maxptr,ounit common/chrbuf/ p,pmax,pmin,outptr,maxptr,ounit c--- functions integer wdthof c--- local varables integer chl,cho,chh,chi,chcoma,chbl,opos,nwid,lhmax,i c--- needed characters data chh,chi,chl,cho,chcoma,chbl /8,9,12,15,45,37/ opos = slbrk - 3 if(hiend) goto 10 call putchr (opos,chl,err) call putchr (0,cho,err) goto 20 10 call print call putchr (opos,chh,err) call putchr (0,chi,err) 20 call putchr (slbrk,chbl,err) if(err .ne. 0) goto 999 nwid = max0 (wdthof(iw(from)),wdthof(iw(to)) ) lhmax = pmax - nwid - 2 do 40 i = from, to call putnum (0,iw(i),nwid,err) call putchr (0,chcoma,err) call putchr (0,chbl,err) if(outptr .lt. lhmax) goto 30 call print call putchr (slbrk,chbl,err) 30 if(err .ne. 0) goto 999 40 continue c--- but dont print the final comma opos = maxptr - 1 call putchr (opos,chbl,err) call print if(.not. hiend) call print 999 return end subroutine depthp (w,iw,n,pt1,pt2,cut,iadjh,hi,rank, + medyet,slbrk,err) c--- compute and print the depth for the current line logical medyet integer n,pt1,pt2,cut,iadjh,hi,rank,slbrk,err integer iw(n) real w(n) c--- W(): holds the n sorted data values c IW(): holds the scaled version of W() c pt1,pt2: pointers into the IW() amd W(). On entry c pt1 = pt2 point to the first data value not yet printed. c On exit, pt1 points to the first data value on the next line, c pt2 is unchanged. c CUT: the largest value on the current (pos) line, or the smallest c value above the current (neg) line c IADJH: points to the high adjacent value in W() and IW() c HI: is the greatest value being displayed c RANK: a running total of the rank from the low end. On exit, c rank is updated to include the count for the current line. c MEDYET: is a logical flag, set .true. when the median value has c been processed. c SLBRK: is the character position on the page of the blank column c between the stem and the leaves. c--- functions integer intfn,wdthof c--- local variables integer chlpar,chrpar,lefcnt,ptz,depth,nwid,opos,ptx integer p(130),pmax,pmin,outptr,maxptr,ounit common/chrbuf/p,pmax,pmin,outptr,maxptr,ounit data chlpar,chrpar/43,44/ ptx = pt1 do 90 pt1=ptx,iadjh if(iw(pt1) .gt. cut) goto 110 if((cut .ge. 0) .and. (iw(pt1) .eq. cut)) goto 110 90 continue c--- last data value if we fall thru here -- point past it for consistency 100 pt1 = iadjh + 1 goto 140 110 if(cut .ne. 0) goto 140 c--- zero cut: if all .le. 0, all zeroes go on "-0" stem if(hi .le. 0) goto 100 c--- both +0 and -0 stems -- share the zeroes between them c first check for numbers rounded to zero -- true -0s do 115 ptz=pt1,n if(w(ptz) .ge. 0.0) goto 117 115 continue 117 pt1 = ptz do 120 ptz=pt1,n if(w(ptz) .gt. 0.0) goto 130 120 continue 130 pt1 = pt1 + intfn(float(ptz - pt1)/2.0,err) c--- compute and print depth 140 lefcnt = pt1 - pt2 rank = rank +lefcnt c--- where is the median ? if(.not. medyet) goto 150 c--- case 1: past the median depth = n - (rank - lefcnt) goto 180 150 if(float(rank) .ne. float(n)/2.0) goto 160 c--- case 2: median falls between stems at this point medyet = .true. goto 170 160 if(float(rank) .lt. float(n+1)/2.0) goto 170 c--- case 3: median is on the current line nwid = wdthof (lefcnt) opos = slbrk - 7 - nwid call putchr (opos,chlpar,err) call putnum (0,lefcnt,nwid,err) call putchr (0,chrpar,err) medyet = .true. goto 999 c--- case 4: not up to median yet 170 depth = rank c--- print the depth, if it hasnt been done yet 180 nwid = wdthof (depth) opos = slbrk - 6 - nwid call putnum (opos,depth,nwid,err) 999 return end subroutine stemp (stem,linwid,negnow,slbrk,err) c--- compute and print the stem logical negnow integer stem,linwid,slbrk,err c--- on entry: c STEM is the inner (near zero) edge of the current line c LINWID is the number of possible different leaf digits c NEGNOW is .true. if the current line is negative c SLBRK is the character position on the page of the blank c column between stems and leaves integer p(130),pmax,pmin,outptr,maxptr,ounit common/chrbuf/p,pmax,pmin,outptr,maxptr,ounit c--- functions integer wdthof c--- local variables integer ch0,chbl,chplus,chmin,chstar,chpt integer nstem,lefdig,nwid,opos,ochr,i,ch5stm(5) data ch0/27/ data chbl,chplus,chmin,chstar,chpt/37,39,40,41,46/ data ch5stm(1),ch5stm(2),ch5stm(3),ch5stm(4)/41,20,6,19/ data ch5stm(5)/46/ nstem = stem/10 lefdig = iabs(stem - nstem*10) nwid = wdthof(nstem) c--- case: how many possible digits/line (= LINWID) if(linwid .ne. 2) goto 260 c--- case1: 2 possible digits/line; 5 lines/stem if(nstem .ne. 0) goto 200 opos = slbrk - 4 if(negnow) call putchr (opos,chmin,err) if(.not. negnow) call putchr(opos,chplus,err) opos = opos + 1 goto 210 200 opos = slbrk - nwid - 2 210 call putnum (opos,nstem,nwid,err) i = lefdig/2 + 1 ochr = ch5stm(i) call putchr (0,ochr,err) goto 990 260 if(linwid .ne. 5) goto 290 c--- case 2: 5 possible digits/line; 2 line/stem opos = slbrk - nwid - 1 if(nstem .ne. 0) goto 270 opos = slbrk - 3 if(negnow) call putchr (opos,chmin,err) if(.not. negnow) call putchr( opos,chplus,err) 270 opos = slbrk - nwid - 1 call putnum (opos,nstem,nwid,err) if(lefdig .lt. 5) call putchr (0,chstar,err) if(lefdig .ge. 5) call putchr (0,chpt,err) goto 990 c--- case 3: 10 possible digits/leaf 1 line/stem 290 if(linwid .eq. 10) goto 300 c--- illegal value -- nice numbers bad ? err = 12 goto 999 300 if(( nstem .ne. 0) .or. .not. negnow) goto 310 opos = slbrk - 3 call putchr (opos,chmin,err) call putchr (0,ch0,err) goto 990 310 opos = slbrk - nwid - 1 call putnum (opos,nstem,nwid,err) 990 call putchr (slbrk,chbl,err) 999 return end subroutine sltitl (unit,err) c--- print the title for a stem-and-leaf display integer err real unit c--- ON ENTRY: unit is the leaf digit unit c--- note that this routine can be modified to print the name of the c batch being displayed if that name is known common/chario/ chars,cmax, + cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk, + chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv, + chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6, + ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh, + chlpar,chrpar,chcoma,chpt integer p(130),pmax,pmin,outptr,maxptr,ounit common/chrbuf/p,pmax,pmin,outptr,maxptr,ounit integer chars(46),cmax integer cha,chb,chc,chd,che,chf,chg,chh,chi,chj,chk, + chl,chm,chn,cho,chp,chq,chr,chs,cht,chu,chv, + chw,chx,chy,chz,ch0,ch1,ch2,ch3,ch4,ch5,ch6, + ch7,ch8,ch9,chbl,cheq,chplus,chmin,chstar,chslsh, + chlpar,chrpar,chcoma,chpt c--- functions integer intfn,wdthof c--- local variables integer iexpt,owid,num,i write(ounit,5000) unit 5000 format( ' STEM AND LEAF DISPLAY ----- LEAF DIGIT UNIT = ',F9.4) call putchr (0,chbl,err) call putchr (0,chbl,err) call putchr (0,ch1,err) call putchr (0,chbl,err) call putchr (0,chbl,err) call putchr (0,ch2,err) call putchr (0,chbl,err) call putchr (0,chbl,err) call putchr (0,chr,err) call putchr (0,che,err) call putchr (0,chp,err) call putchr (0,chr,err) call putchr (0,che,err) call putchr (0,chs,err) call putchr (0,che,err) call putchr (0,chn,err) call putchr (0,cht,err) call putchr (0,chs,err) call putchr (0,chbl,err) iexpt = intfn(alog10(unit),err) if(iexpt .ge. 0) goto 200 if(iexpt .eq. (-1)) goto 100 c--- unit .le. 0.01 iexpt = iabs(iexpt) - 2 call putchr (0,ch0,err) call putchr (0,chpt,err) if(iexpt .eq. 0) goto 30 do 20 i=1,iexpt call putchr (0,ch0,err) 20 continue 30 call putchr (0,ch1,err) call putchr (0,ch2,err) goto 900 100 call putchr (0,ch1,err) call putchr (0,chpt,err) call putchr (0,ch2,err) goto 900 c--- unit .ge. 1.0 200 num = 12 * intfn(unit,err) owid = wdthof(num) call putnum (0,num,owid,err) call putchr (0,chpt,err) c--- wrap up 900 if(err .ne. 0) goto 999 call print write (ounit,5010) 5010 format(/) 999 return end subroutine rgcomp (x,y,l,mu,sigma,yhat,drr,mhat,shat,err) integer l,err real x(l),y(l),yhat(l),drr(l),mu,sigma,mhat,shat c--- Perform the computations for a suspended rootogram. c X(I) ... X(L) are the bin boundaries and Y(1) ... Y(L) c are the bin counts. The count Y(I) corresponds to the bin whose c right boundary is X(I). The bin whose right boundary is X(1) c is open to the left. Also X(l) is now used, so that Y(L) counts c all data values to the right of X(L-1). c A Gaussian comparison curve is used, and its center and scaled c are determined by the hinges of the data (found by linear interpolation). c If sigma is not equal to zero, then the fitting process is skipped c and the values of MU and SIGMA passed in are used for the comparison c curve. If SIGMA is equal to zero, the values of both MU and SIGMA c are ignored. c On exit, MHAT contains the fitted mean of the gaussian comparison c curve, and SHAT contains the fitted standard deviation, YHAT() contains c the L fitted counts, and DRR() contains the double-root residuals. integer i,k,lp1,lp1mi real d,hl,hu,p,pl,t,tn,yh if(l .ge. 3) goto 5 err = 91 goto 999 5 k = l - 1 tn = 0.0 do 10 i =1,l tn = tn + y(i) 10 continue c--- if mu and sigma were specified, dont bother to fit them from c the data. Cue is non-zero SIGMA. if(sigma .gt. 0.0) goto 80 d = 0.5 *(1.0 + aint(0.5*(tn+1.0))) c--- if lower hinge falls in left-open bin, error. if(d .gt. y(1)) goto 20 err = 92 goto 999 20 t = y(1) do 30 i=2,k t = t + y(i) if(t .ge. d) goto 40 30 continue c--- lower hinge falls in right-open bin -- error. err =92 goto 999 c--- find lower hinge by interpolation 40 t = t - y(i) hl = x(i-1) + (x(i) - x(i-1)) * (d - t - 0.5)/y(i) c--- do same thing for upper hinge if(d .gt. y(l)) goto 50 err = 92 goto 999 50 t = y(l) lp1 = l + 1 do 60 i =2,k lp1mi = lp1 -i t = t + y(lp1mi) if(t .ge. d) goto 70 60 continue err = 92 goto 999 70 t = t - y(lp1mi) hu = x(lp1mi) - (x(lp1mi) - x(lp1mi-1)) * (d - t - 0.5)/y(lp1mi) c--- use mhat = mid-hinge for centering and shat = (h-spread)/1.349 c for scale. mhat = (hl + hu)/2.0 shat = (hu - hl)/1.349 goto 90 c--- skip to here if mu and sigma were specified 80 mhat = mu shat = sigma 90 pl = 0.0 do 100 i=1,k p = gau((x(i) - mhat)/shat) yh = tn * (p -pl) yhat(i) = yh drr(i) = sqrt(2.0+4.0*y(i)) - sqrt(1.0+4.0*yh) if(y(i) .eq. 0.0) drr(i) = 1.0 - sqrt(1.0+4.0*yh) pl = p 100 continue yh = tn * (1.0 - pl) yhat(l) = yh drr(l) = sqrt(2.0+4.0*y(l)) - sqrt(1.0+4.0*yh) if(y(l) .eq. 0.0) drr(l) = 1.0 - sqrt(1.0+4.0*yh) 999 return end subroutine rgprnt (y,l,yhat,drr,err) integer l,err real y(l),yhat(l),drr(l) c--- Print, bin by bin, the observed count, the raw residual, the c double root residual and an abbreviated display of the double c root residual. Y(1) ... Y(L) contains the fitted counts, and DRR c contains the double root residuals. integer bl,bo,dot,i,j,min,nbl,nmin,npl,pl,star real res integer floor integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit data bl,dot,min,pl,star /1h ,1h.,1h-,1h+,1h*/ c--- is print line wide enought to hold the numbers if(pmax .ge. 30) goto 10 err = 93 goto 999 c--- print line may be adequate for numbers but not display 10 if(pmax .ge. 65) goto 30 err = 94 c--- print only the observed counts and the two types of residuals write (ounit,5010) 5010 format (1x,'BIN',3x,'COUNT',3x,'RAWRES',3x,'DRRES'/) do 20 i=1,l res = y(i) - yhat(i) write (ounit,5020) i,y(i),res,drr(i) 20 continue 5020 format(1x,i3,2x,f6.1,4x,f5.1,4x,f5.2) goto 999 c--- print the table and the display 30 write (ounit,5030) 5030 format(1x,'BIN',3x,'COUNT',3x,'RAWRES',4x,'DRRES',7x, + 'SUSPENDED ROOTOGRAM'/) do 120 i=1,l res = y(i) - yhat(i) if (drr(i) .ne. 0.0) goto 40 write(ounit,5040) i,y(i),res,drr(i) 5040 format(1x,i3,2x,f6.1,4x,f5.1,4x,f5.2,8x,'.',20x,'.') goto 120 40 if(drr(i) .gt. 0.0) goto 80 c--- handle lines with negative drr -- there are four cases c -s end in * to indicate overflow c -s overwrite dot but fit on line c no blanks between dot and -s, and c at least one blank between dot and -s nmin = -floor(5.0*drr(i)) if(nmin .gt. 10) goto 60 if(nmin .lt. 10) goto 50 write (ounit,5050) i,y(i),res,drr(i), + (bl,j=1,5),dot,(min,j=1,10), + (bl,j=1,10), dot 5050 format(1x,i3,2x,f6.1,4x,f5.1,4x,f5.2,3x,32a1) goto 120 50 nbl = 10 - nmin write (ounit, 5050) i,y(i),res,drr(i), + (bl,j=1,5),dot,(bl,j=1,nbl),(min,j=1,nmin), + (bl, j=1,10),dot goto 120 60 bo = bl if(nmin .le. 15) goto 70 nmin = 15 bo = star 70 nbl = 16 - nmin write (ounit,5050) i,y(i),res,drr(i), + (bo,j=1,nbl),(min,j=1,nmin),(bl,j=1,10),dot goto 120 c--- handle lines with positive drr -- there are four cases c +s end in * to indicated overflow c +s overwrite dot but fit on line c no blanks between dot and +s and c at least 1 blank between dot and +s 80 npl = -floor(-5.0*drr(i)) if(npl .gt. 10) goto 100 if(npl .lt. 10) goto 90 write (ounit,5050) i,y(i),res,drr(i), + (bl,j=1,5),dot,(bl,j=1,10),(pl,j=1,10),dot goto 120 90 nbl = 10 - npl write (ounit,5050) i,y(i),res,drr(i), + (bl,j=1,5),dot,(bl,j=1,10),(pl,j=1,npl),(bl,j=1,nbl),dot goto 120 100 bo = bl if(npl .le. 15) goto 110 npl = 15 bo = star 110 nbl = 16 - npl write (ounit,5050) i,y(i),res,drr(i), + (bl,j=1,5),dot,(bl,j=1,10),(pl,j=1,npl),(bo,j=1,nbl) 120 continue write(ounit,5060) 5060 format(/1x,' IN DISPLAY, VALUE OF ONE CHARACTER IS .2', + 7x,'00'/) 999 return end subroutine plot (y,x,n,wy,wx,linset,chrset,xmin,xmax, + ymin,ymax,err) c--- Plot the N ordered pairs (X(I),Y(I)) using a condensed plot. c Condensed plotting used the plotting symbol to indicate c the fine detail of numberical spacing. As a result more precision c can be conveyed in fewer lines. Multiple points falling at the same c plot position are not indicated, however -- the most extreme (in Y) c point will be selected for display. c X() and Y() are not modified by the program. Work is done using c the work arrays WX() and WY() subblied fy the calling program. c The details of plot format are determined by the parameters in the c calling sequence. LINSET specifies the maximum number of lines to c be used. CHRSET specifies how many different codes can be used on c each line. If either of these is zero, the program defaults to c 10 lines and 10 character codes (0 thru 9). c XMIN and XMAX specify the range of the x-values to be plotted. c YMIN and YMAX specify the range of the y-values to be plotted. c For either pair, if they are set equal by the calling program, c the program defaults to using the adjacent values in each dimension. c This option is almost always preferred for exploratory plots. integer n,linset,chrset,err real y(n),x(n),wy(n),wx(n),xmin,xmax,ymin,ymax integer p(130),pmax,pmin,outptr,maxptr,ounit common/chrbuf/ p,pmax,pmin,outptr,maxptr,ounit c--- functions integer intfn,floor,wdthof c--- call subroutines nposw,print,qsort2,putchr,putnum,yinfo c--- loacal variables integer chl,chm,chp,chr,ch0,ch9,chplus,chmin integer chrpar,chstar integer lines,chrs,maxl,xposns,iadjl,iadjh,nn,lftpsn integer lnsfrz,lineno,ptr,nwid,proom,ochar,opos,ychar integer oposx,lnflor,label,i real hh,hl,med,step,top,bottom,left,right real adjxl,ahjxh,adjyl,adjyh,xfract,xunit,xnpw,yfract real yunit,ynpw,ylabel,xval,syval,nionos(9) logical negnow,markzs data chl,chm,chp,chr,ch0,ch9/12,13,16,18,27,36/ data chplus,chmin,chstar,chrpar/39,40,41,44/ data nn,nionos(1),nionos(2),nionos(3) /9,1.0,1.5,2.0/ data nionos(4),nionos(5),nionos(6) /2.5,3.0,4.0/ data nionos(7),nionos(8),nionos(9) /5.0, 7.0, 10.0/ data markzs /.false./ if(n .ge. 5) goto 10 err = 41 goto 999 10 lftpsn = pmin + 6 lines = 10 chrs = 10 if(linset .eq. 0 .or. chrset .eq. 0) goto 30 lines = linset chrs = chrset err = 42 if(lines .lt. 5 .or. lines .gt. 40) goto 999 if(chrs .lt. 1 .or. chrs .gt. 10) goto 999 err = 0 c--- set up scales and plot boundary information 30 lftpsn = pmin + 6 proom = pmax - lftpsn + 1 do 40 i=1,n wx(i) = x(i) 40 continue if(xmin .ge. xmax) goto 45 call yinfo (wx,n,med,hl,hh,adjxl,adjxh,iadjl,iadjh,step,err) if(err .ne. 0) goto 999 if(adjxl .lt. adjxh) goto 50 c--- if x-adjacent values equal, try using the extremes 45 adjxl = wx(1) adjxh = wx(n) err = 44 if(adjxl .ge. adjxh) goto 999 err = 0 50 call nposw (adjxh,adjxl,nionos,nn,proom,.false.,xposns, + xfract,xunit,xnpw,err) c--- scale y --- sort(x,y) paired on y do 60 i=1,n wx(i) = x(i) wy(i) = y(i) 60 continue call qsort2(n,wy,wx) if(ymin .ge. ymax) goto 65 call yinfo (wy,n,med,hl,hh,adjyl,adjyh,iadjl,iadjh,step,err) if(err .ne. 0) goto 999 goto 68 65 adjyl = wy(1) adjyh = wy(n) err = 45 if(adjyl .ge. adjyh) goto 999 err = 0 68 maxl = lines call nposw (adjyh,adjyl,nionos,nn,maxl,.true.,lines, + yfract,yunit,ynpw,err) if(err .ne. 0) goto 999 if(yfract .ne. 10.0) goto 70 yfract = 1.0 yunit = yunit*10.0 c--- find nice plot edges -- round away from center of plot 70 lnsfrz = -floor(-adjyh/ynpw) top = float(lnsfrz)*ynpw bottom = float(floor(adjyl/ynpw))*ynpw left = float(floor(adjxl/xnpw))*xnpw right = float(-floor(-adjxh/xnpw))*xnpw c--- print scrawl write (ounit,9070) lines,chrs 9070 format (15x,i3,7H LINE, , i3,15H CHARACTER PLOT) write (ounit,9080) bottom,top,ynpw,left,right,xnpw 9080 format (15x, 8H Y FROM , f12.6, 4H TO , f12.6, 7H STEP , + f12.6/,15x, 8H X FROM , f12.6, 4H TO , f12.6, + 7H STEP , F12.6//) c--- initialize for plotting -- one line too high c LNSFRZ counts # lines away from zero == +00 and -00 are 0 lines away ylabel = float(lnsfrz) * yfract lnflor = lnsfrz negnow = .false. if(top .gt. 0.0) goto 80 lnsfrz = lnsfrz + 1 negnow = .true. 80 lineno = 0 ptr = n+1 c--- start a new line of the plot 90 lnflor = lnflor - 1 lineno = lineno + 1 opos = pmin if(lnsfrz .gt. 0 .or. negnow) goto 95 c--- just went negative negnow = .true. goto 97 95 lnsfrz = lnsfrz - 1 ylabel = ylabel - yfract 97 continue c--- print the line label if(.not.negnow) call putchr(opos,chplus,err) if(negnow) call putchr(opos,chmin,err) if(ylabel .ne. 0.0) goto 120 opos = pmin + 3 call putchr (opos,ch0,err) call putchr (0,ch0,err) if((chrs .gt. 1) .and. ((lines-lineno) .ge. 3)) markzs = .true. opos = pmin + 5 call putchr (opos,chrpar,err) if(.not. markzs ) goto 111 do 100 i=1,5 if(.not. negnow) call putchr(0,chmin,err) if(negnow) call putchr(0,chmin,err) 100 continue oposx = pmax - 5 do 110 opos = oposx,pmax if(.not. negnow) call putchr (0,chmin,err) if( negnow) call putchr (0,chmin,err) 110 continue 111 continue c--- print non-zero label 120 label = intfn (10.0 * abs(ylabel),err) if(err .ne. 0) goto 999 nwid = wdthof (label) opos = pmin + 5 - nwid call putnum (opos,label,nwid,err) if(err .ne. 0) goto 999 c--- get next data point 125 ptr = ptr - 1 if(ptr .le. 0) goto 135 xval = wx(ptr) syval = wy(ptr)/ynpw if(intfn(syval,err) .gt. lnflor) goto 140 if(intfn(syval,err) .eq. lnflor .and. syval .ge. 0.0) goto 140 c--- time to start next line c if this is the last line, print it anyway and use "m" for low no. 130 if(lineno .eq. lines) goto 140 c--- back up the pointer ptr = ptr + 1 c--- wrap up line 135 if(err .ne. 0) goto 999 call print c--- and start a new line goto 90 c--- get y-character 140 ychar = ifix(abs(syval - float(lnsfrz))*float(chrs)) ochar = ch0 + ychar if(ochar .ge. ch0 .and. ochar .le. ch9) goto 145 ochar = chp if(lineno .eq. lines) ochar = chm c--- get x-position 145 opos = pmin + 5 + intfn((xval - left)/xnpw,err) + 1 if(xval .ge. left) goto 150 opos = pmin + 6 if(ochar .lt. ch0 .or. ochar .gt. ch9) goto 147 ochar = chl goto 160 147 ochar = chstar goto 160 150 if(xval .le. right) goto 160 opos = pmax if(ochar .lt. ch0 .or. ochar .gt. ch9) goto 157 ochar = chr goto 160 157 ochar = chstar 160 continue call putchr (opos,ochar,err) if(err .ne. 0) goto 999 if(ptr .gt. 1) goto 125 call print 999 return end subroutine boxes (y,n,gsub,ng,line3,notch,sorty,err) c--- print adjacent boxplots on a single scale for all variables in y(). integer n,ng,err integer gsub(n) real y(n),sorty(n) logical line3,notch c--- Y() contains data. GSUB() contains integers between 1 and ng c identifying the data set each element of Y() belongs to. c This data structure is consistent with the sparse matrix c format used for storing matrices in other programs. The use of c the vector GSUB() is meant to suggest boxplots of either the c rows or the columns a matrix stored in this manner. c If LINE3 is .true. all boxplots will be full 3-line boxplots. c If LINE3 is .false., one-line boxplots will be printed. c Scaling of these plots is to the extremes of the entire combined c data batch. The details of each box, including outlier identification, c are determined for each batch individually. integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit integer nn,npmax,npos,lpmin,spmin integer chrpar, lblw, opos,i,j,k real nionos(9), fract,unit,npw,lo,hi integer wdthof c--- call subroutines boxp,nposw,putchr,putnum data nn,nionos(1),nionos(2),nionos(3)/9,1.0,1.5,2.0/ data nionos(4),nionos(5),nionos(6)/2.5,3.0,4.0/ data nionos(7),nionos(8),nionos(9)/5.0,7.0,10.0/ data chrpar/44/ c--- check for at least two data values. Otherwise highest and lowest c will be equal and scaling will fail. if(n .gt. 1) goto 5 err = 31 goto 999 5 lpmin = pmin + 7 lo = y(1) hi = y(n) do 10 i =1,n if(lo .gt. y(i)) lo = y(i) if(hi .lt. y(i)) hi = y(i) 10 continue c--- scale to the extremes npmax = pmax - lpmin + 1 call nposw (hi,lo,nionos,nn,npmax,.false.,npos,fract, + unit,npw,err) if(err .ne. 0) goto 999 c--- now print all the boxes. Data sets identified by their codes c in GSUB() if(ng .gt. 1) goto 17 do 15 k=1,n sorty(k) = y(k) 15 continue call boxp (sorty,n,line3,notch,lo,hi,npw,err) goto 999 17 spmin = pmin do 30 i=1,ng k=0 do 20 j=1,n if(gsub(j) .ne. i) goto 20 k = k+1 sorty(k) = y(j) 20 continue pmin = spmin lblw = wdthof(i) opos = pmin + 5 - lblw call putnum (opos,i,lblw,err) opos = pmin + 6 call putchr (opos,chrpar,err) if(err .ne. 0) goto 999 pmin = lpmin call boxp (sorty,k,line3,notch,lo,hi,npw,err) if(err .ne. 0) goto 999 30 continue pmin = spmin 999 return end subroutine boxp (sorty,n,line3,notch,lo,hi,npw,err) c--- print a boxplot of the data in SORTY() integer n,err real sorty(n),lo,hi,npw logical line3,notch c--- Plot scaling has been done by the calling program with needed c information passed in as lo (the low extreme), hi (the high extreme) c and npw (the nice position width for plotting). c Typically this will be one of several boxplots scaled and printed c together. c If LINE3 is .true. a 3-line boxplot (full boxes) is printed. c If not, the simple one-line boxplot is printed. Both convey the c same information, but the 3-line version may look nicer. c If NOTCH is .true. a confidence interval around the median is indicated c with parenthesis integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit integer pltpos c--- calls subroutine boxtop,print,putchr,yinfo integer i,iadjl,iadjh,ifrom,ito,lpmax,lpmin integer opos,chi,cho,chstar,chmin,chplus,chrpar,chlpar real med,hl,hh,adjl,adjh,step real floatn,nstep,lnotch,hnotch,ofencl,ofench data chi,cho,chplus,chmin,chstar /9,15,39,40,41/ data chlpar,chrpar /43,44/ lpmax = pmax lpmin = pmin call yinfo (sorty,n,med,hl,hh,adjl,adjh,iadjl,iadjh,step,err) if(err .ne. 0) goto 999 floatn = float(n) nstep = 1.7 * (1.25*(HH- HL)/(1.35*sqrt(floatn))) lnotch = med - nstep hnotch = med + nstep c--- print top of box, if 3-line version if(line3) call boxtop (lo,hi,hl,hh,npw,err) if(err .ne. 0) goto 999 c--- fill center line of display -- note careful heirarchy c of over printing. Last placed character is only one to appear c--- mark whiskers ifrom = pltpos (adjl,lo,hi,npw,err) ito = pltpos (hl,lo,hi,npw,err) - 1 if(ifrom .gt. ito) goto 21 do 20 i=ifrom,ito call putchr (i,chmin,err) 20 continue 21 continue ifrom = pltpos (hh,lo,hi,npw,err) ito = pltpos (adjh,lo,hi,npw,err) if(ifrom .gt. ito) goto 31 do 30 i=ifrom,ito call putchr (i,chmin,err) 30 continue 31 continue c--- mark low outliers, if any if(iadjl .eq. 1) goto 41 ofencl = hl - 2.0*step ito = iadjl - 1 do 40 i=1,ito opos = pltpos (sorty(i),lo,hi,npw,err) if(sorty(i) .lt. ofencl) call putchr (opos,cho,err) if(sorty(i) .ge. ofencl) call putchr (opos,chstar,err) 40 continue 41 continue c--- mark high outliers, if any if(iadjh .eq. n) goto 51 ofench = hh + 2.0*step ifrom = iadjh + 1 do 50 i=ifrom,n opos = pltpos (sorty(i),lo,hi,npw,err) if(sorty(i) .gt. ofench) call putchr (opos,cho,err) if(sorty(i) .le. ofench) call putchr (opos,chstar,err) 50 continue 51 continue c--- mark hinges, notches, and median opos = pltpos (hl,lo,hi,npw,err) call putchr (opos,chi,err) opos = pltpos (hh,lo,hi,npw,err) call putchr (opos,chi,err) opos = pltpos (lnotch,lo,hi,npw,err) if(notch) call putchr (opos,chlpar,err) opos = pltpos (hnotch,lo,hi,npw,err) if(notch) call putchr (opos,chrpar,err) opos = pltpos (med,lo,hi,npw,err) call putchr (opos,chplus,err) c--- print the boxplot if(err .ne. 0) goto 999 call print c--- print the bottom of the box if(line3) call boxtop (lo,hi,hl,hh,npw,err) 999 return end subroutine boxtop (lo,hi,hl,hh,npw,err) real lo,hi,hl,hh,npw integer err c--- print the top or bottom of a boxplot display c--- hi and lo are adges of the plotting region used by the pltpos c function. c hl and hh are the low and high hinges c mpw is the nice position width set by the plot scaling routines integer i,ifrom,ito,chmin integer pltpos data chmin/40/ ifrom = pltpos (hl,lo,hi,npw,err) ito = pltpos (hh,lo,hi,npw,err) if(ifrom .gt. ito) goto 11 do 10 i=ifrom,ito call putchr (i,chmin,err) 10 continue 11 continue if(err .eq. 0) call print 999 return end integer function pltpos (x,lo,hi,npw,err) c--- find the position corresponding to c on plot bounded c between lo and hi and scaled according to npw real x,lo,hi,npw integer err integer intfn integer p(130),pmax,pmin,outptr,maxptr,ounit common /chrbuf/ p,pmax,pmin,outptr,maxptr,ounit pltpos = intfn ((x-lo)/npw,err) + pmin if(pltpos .lt. pmin) pltpos = pmin if(pltpos .gt. pmax) pltpos = pmax return end subroutine rline (x,y,n,resid,work,nsteps,slope,level,lls, c lus,trace,lhslop,rhslop,hsrtio,err) integer n,nsteps,err real x(n),y(n),resid(n),work(n),slope,level,lls,lus real lhslop,rhslop,hsrtio logical trace c--- For the data (x(1),y(1)), ... , (x(n),y(n)), fit the straight line c y = level + slope*x + resid c by the "resistant line" technique. Iterates for nsteps steps or c until the slope is correct to four digits. 1/tol specifies the c number of digits required. If convergence is not attained after c nsteps steps, llu and lls will return the last lower and upper c bounds on the correct slope -- otherwise they will return zero. c This method will not work for n .le. 5 and it will not be fully c resistant for n .le. 7. If several x-values are tied, n should c be still larger to guarantee resistance. The program also computes c the approxiamte slope of the left half and the right half of the c data in lhslop and rhslop. Their ratio, returned in hsrtio, is a c measure of the straightness of the x-y relationship. c If trace is .true. on entry, the halfslope entry will be printed and c a report will be printed after each step of the iteration. c--- common common /numbrs/ epsi,maxint real epsi,maxint common /chrbuf/ p,pmax,omin,outptr,maxptr,ounit integer p(130),pmax,pmin,outptr,matptr,ounit c--- local variables integer i,mpt1,mpt2,mpt3,n1,n2,n3,l3rd,r3rd integer mxlo3,mnhi3,from,to,stepno,mptx,mpty real x1,x2,x3,y1,y2,y3,xmed,dslope,tol,tol1 real slope1,slope2,slope3,deltx,dr1,dr2,dr3 real oldds,ddr,numtor,denom,dsd2 c--- functions real rl3med, deltr, median c--- 1/tol specifies number of reliable slope digits required tol = 1.0e-4 lls = 0.0 lus = 0.0 if (n.gt.5) goto 5 err = 51 goto 999 5 if (nsteps.gt.0) goto 10 err = 52 goto 999 c--- divide into thids on x c first check for ties 10 call qsort2 (n,x,y) err = 0 if(err.ne.0) goto 999 mpt2 = (n/2) + 1 mpt1 = n - mpt2 + 1 xmed = (x(mpt1) + x(mpt2))/2.0 c--- look for first value not tied with median. it is the maximum c possible for low third cut mxlo3 = mpt2 20 mxlo3 = mxlo3 - 1 if(x(mxlo3) .ne. xmed) goto 30 if(mxlo3 .gt. 1) goto 20 c--- fall though here is all tied from low end to median mxlo3 = 0 c--- look for minumum possible high third cut 30 mnhi3 = mpt1 40 mnhi3 = mnhi3 + 1 if(x(mnhi3) .ne. xmed) goto 60 if(mnhi3 .lt. n) goto 40 c--- fall thorugh here if all tied from median to high end err = 53 goto 999 c--- only two "thirds" 50 mnhi3 = mxlo3 +1 goto 70 60 if(mxlo3 .ne. 0) goto 70 c--- low third empty mxlo3 = mnhi3 - 1 70 continue c--- now place the thirds c get favored low end split mpt1 = (n+1)/3 x1 = x(mpt1) c--- dont split ties c favor larger outer thirds l3rd = mxlo3 if(mpt1 .gt. mxlo3) goto 90 l3rd = mpt1 80 l3rd = l3rd +1 if(x(l3rd) .eq. x1) goto 80 l3rd = l3rd - 1 c--- now the high third 90 mpt3 = n - mpt1 + 1 x3 = x(mpt3) c--- dont split ties c favor larger outer thirds r3rd = mnhi3 if(mpt3 .le. mnhi3) goto 110 r3rd = mpt3 100 r3rd = r3rd -1 if(x(r3rd) .eq. x3) goto 100 r3rd = r3rd +1 110 continue c--- now l3rd and r3rd point to inner edges of outer thirds c--- check if thirds are big enough for resistance n1 = l3rd n3 = n - r3rd + 1 n2 = n - n1 - n3 if ((n1. gt. 2) .or. (n3 .gt. 2)) goto 120 c--- if n = 7 and split is 2 - 3 - 2, stick with it if ((n1 .eq. 2) .and. (n2 .eq. 3) .and. (n3.eq.2)) goto 140 err = 54 goto 999 120 if ((n1 .gt. 2) .and. (n3.gt. 2)) goto 140 c--- only 2 thirds are big enough -- regroup and work with 2 if(n1 .le. 2) l3rd = r3rd -1 if(n3 .le. 2) r3rd = l3rd + 1 130 n1 = l3rd n2 = 0 n3 = n - r3rd + 1 x2 = 0.0 y2 = 0.0 140 continue c--- set up for fitting c--- get x medians mpt1 = (n1+1)/2 mpt2 = (n2+1)/2 mpt3 = (n3+1)/2 mpty = n1 - mpt1 + 1 x1 = (x(mpt1) + x(mpty))/2.0 mptx = n1 + mpt2 mpty = n1 + n2 - mpt2 + 1 if(n2 .ne. 0) x2 = (x(mptx) + x(mpty))/2.0 mptx = n1 + n2 + mpt3 mpty = n - mpt3 + 1 x3 = (x(mptx) + x(mpty))/2.0 deltx = x3 - x1 if(abs(deltx) .lt. epsi) deltx = sign(epsi,deltx) c--- y-medians y1 = rl3med (y,n,1,l3rd,work,err) from = l3rd + 1 to = r3rd - 1 if (n2 .ne. 0) y2 = rl3med (y,n,from,to,work,err) y3 = rl3med (y,n,r3rd,n,work,err) if (err .ne. 0) goto 999 c--- compute half-slope ratio to check straightness of y on x c report if trace is .true. else just return results if( n2 .eq. 0) goto 170 lhslop = (y2 - y1)/(x2 - x1) rhslop = (y3 - y2)/(x3 - x2) if(abs(lhslop) .gt. epsi) goto 160 hsrtio = 0.0 goto 170 160 hsrtio = rhslop/lhslop if(trace) write(ounit,5002) lhslop,rhslop,hsrtio 5002 format(1x,'straightness check'/1x,'left half-slope =', c f12.6,' right half-slope = ',f12.6/10x,'ratio =',f12.6//) 170 continue c--- first 2 slopes without iterating stepno = 1 slope1 = (y3 - y1)/deltx dr1 = deltr (x,y,n,resid,l3rd,r3rd,slope1,work,err) if(err.ne.0) goto 999 dslope = dr1/deltx if(trace) write(ounit,5000) stepno,slope1 5000 format(1x,'slope ',i3,': ',f12.6) stepno = 2 slope2 = slope1 + dslope slope3 = slope1 180 dr2 = deltr (x,y,n,reaid,l3rd,r3rd,slope2,work,err) if(err .ne. 0) goto 999 if(dr2 .eq. 0.0) goto 290 c--- find second slope with opposite sign residual difference if(sign(1.0,dr2) .ne. sign(1.0,dr1)) goto 190 slope1 = slope2 dr1 = dr2 slope2 = slope3 + dslope dslope = dslope + dslope goto 180 190 if (trace) write (ounit,5000) stepno,slope2 adr = abs(dr2) c--- iteration is based upon the algorithm ZEROIN 220 slope3 = slope1 dr3 = dr1 dslope = slope2 - slope1 oldds = dslope 230 if( abs(dr3) .ge. abs(dr2)) goto 240 slope1 = slope2 slope2 = slope3 slope3 = slope1 dr1 = dr2 dr2 = dr3 dr3 = dr1 c--- test convergence 240 if( stepno .ge. nsteps) goto 285 tol1 = 2.0 * epsi * abs(slope2) +0.5*tol dsd2 = .5 * (slope3 - slope2) if(abs(dsd2) .le. tol1) goto 290 if(dr2 .eq. 0.0) goto 290 c--- try again ddr = dr2/dr1 numtor = 2.0*dsd2*ddr denom = 1.0 - ddr if(numtor .gt. 0.0) denom = -denom numtor = abs(numtor) if((2.0*numtor).ge.(3.0*dsd2*denom-abs(tol1*denom))) goto 270 if(numtor .ge. abs(0.5*oldds*denom)) goto 270 oldds = dslope dslope = numtor/denom goto 280 c--- bisect 270 dslope = dsd2 oldds = dslope 280 slope1 = slope2 dr1 = dr2 if( abs(dslope) .gt. tol1) slope2 = slope2 + dslope if( abs(dslope) .le. tol1) slope2 = slope2 + sign(tol1,dsd2) stepno = stepno + 1 if(trace) write(ounit,5000) stepno,slope2 if( err .ne.0) goto 999 if((dr2 * (dr3/abs(dr3))) .gt. 0.0) goto 220 goto 230 c--- ran out of steps 285 lls = amin1(slope1,slope3) lus = amax1(slope1,slope3) goto 999 c--- exit 290 slope = slope2 do 300 i=1,n work(i) = y(i) - slope * x(i) 300 continue call qsort4 (n,work) if(err .ne. 0) goto 999 level = median (work,n) do 310 i=1,n resid(i) = y(i) - slope*x(i) - level 310 continue 999 return end real function rl3med (y,n,from,to,work,err) c--- returns the median of the numbers from y(from) to y(to), inclusive integer n,from,to,err real y(n),work(n) integer i,j real median j=0 do 10 i= from,to j = j+1 work(j) = y(i) 10 continue call qsort4 (j,work) if(err .ne. 0) goto 999 rl3med = median (work,j) 999 return end real function deltr (x,y,n,resid,l3rd,r3rd,slope,work,err) c--- returns the difference between the median residuals in the left c and right 3rds of the data for a line with the specified slope integer n,l3rd,r3rd,err real x(n),y(n),resid(n),work(n),slope integer i real rl3med do 10 i=1,n resid(i) = y(i) - slope*x(i) 10 continue deltr = rl3med(resid,n,r3rd,n,work,err) - c rl3med(resid,n,1,l3rd,work,err) return end subroutine rsm (y,n,smooth,rough,versn,err) c--- main program for nonlinear smoothers integer n,versn,err,i real y(n),smooth(n),rough(n) c--- on entry: c y() is a dta sequence of n values c versn specifies the smoother to be used c versn=1 specifies 3rssh, twice c versn=2 specifies 4253h, twice c--- on exit: c smooth() and rough() contain the smooth and rough c resulting from the smoothing operation. Note that c y(i) = smooth(i) + rough(i) c for each i from 1 to n if (n.gt.6) goto 10 err = 61 goto 999 10 do 20 i=1,n smooth(i) = y(i) 20 continue if(versn.eq.1) call s3rssh(smooth,n,err) if(versn.eq.2) call s4253h(smooth,n,err) c--- compute rough from first smoothing do 30 i=1,n rough(i) = y(i) - smooth(i) 30 continue c--- rerough smoothers ("twicing") if(versn.eq.1) call s3rssh (rough,n,err) if(versn.eq.2) call s4253h (rough,n,err) if(err.ne.0) goto 999 do 40 i=1,n smooth(i) = smooth(i) + rough(i) rough(i) = y(i) - smooth(i) 40 continue 999 return end subroutine s3rssh (y,n,err) c--- smooth y() by 3rssh, twice integer n,err real y(n) logical change call s3r (y,n) change = .false. call split (y,n,change) if(.not. change) goto 10 call s3r (y,n) change = .false. call split (y,n,change) if(change) call s3r (y,n) 10 call hann (y,n) 999 return end subroutine s4253h (y,n,err) c--- smooth by 4253h integer n,err real y(n) real endsav,work(5),save(5) integer nw logical change data nw/5/ change = .false. call s4 (y,n,endsav,work,save,nw,err) if(err .eq. 0) call s2 (y,n,endsav) if(err .eq. 0) call s5 (y,n,work,save,nw,err) if(err .eq. 0) call s3 (y,n,change) if(err .eq. 0) call endpts (y,n) if(err .eq. 0) call hann (y,n) 999 return end subroutine s4 (y,n,endsav,work,save,nw,err) c--- smooth by running medians of 4 integer n,nw,err real y(n),endsav,work(nw),save(nw) real endm1,two data two/2.0/ c--- even length medians offset the output sequence to the high end, c since they cannot be symmetric. endsav is left holding y(n) since c there is no other room for it. y(1) is unchanged. endsav = y(n) endm1 = y(n-1) call runmed (y,n,4,work,save,nw,err) y(2) = (y(1) + y(2))/two y(n) = (endm1 + endsav)/two 999 return end subroutine s2 (y,n,endsav) c--- smooth by running medians (means) of 2. c used to recenter results of running medians of 4. c endsav holds the original y(n). integer n real y(n),endsav integer nm1,i real two data two/2.0/ nm1 = n-1 do 10 i=2,nm1 y(i) = (y(i+1)+y(i))/two 10 continue y(n) = endsav 999 return end subroutine s5 (y,n,work,save,nw,err) c--- smooth by running medians of 5 integer n,nw,err real y(n),work(nw),save(nw) logical change real ymed1,ymed2 change = .false. call medof3 (y(1),y(2),y(3),ymed1,change) call medof3 (y(n),y(n-1),y(n-2),ymed2,change) call runmed (y,n,5,work,save,nw,err) y(2) = ymed1 y(n-1) = ymed2 999 return end subroutine hann (y,n) c--- 3-point smooth by movind averages weighted 1/4,1/2,1/4. integer n real y(n) integer i,nm1 real y1,y2,y3 nm1 = n-1 y2 = y(1) y3 = y(2) do 10 i=2,nm1 y1 = y2 y2 = y3 y3 = y(i+1) y(i) = (y1 + y2 + y2 + y3)/4.0 10 continue 999 return end subroutine s3 (y,n,change) c--- compute running median of 3 on y(i) c sets change .true. if any change made integer n real y(n) logical change real y1,y2,y3 integer nm1 y2 = y(1) y3 = y(2) nm1 = n-1 do 10 i=2,nm1 y1 = y2 y2 = y3 y3 = y(i+1) call medof3 (y1,y2,y3,y(i),change) 10 continue 999 return end subroutine s3r (y,n) c--- compute running medians of 3 integer n real y(n) logical change 10 change = .false. call s3 (y,n,change) if(change) goto 10 call endpts (y,n) 999 return end subroutine medof3 (x1,x2,x3,xmed,change) c--- put the median of x1, x2,x2 in xmed and set change .true. c if the median isnt x2 real x1,x2,x3,xmed logical change real y1,y2,y3 y1 = x1 y2 = x2 y3 = x3 xmed = y2 if ((y2-y1)*(y3-y2) .ge. 0.0) goto 999 change = .true. xmed = y1 if ((y3-y1)*(y3-y2) .gt. 0.0) goto 999 xmed = y3 999 return end subroutine endpts (y,n) c--- estimate smoothed values for both endpoints of the sequence in y() c using the end point extrapolation rule. c all the values in y(0 except the end points have been smoothed. integer n real y(n) real y0,ymed logical change change = .false. c--- left end y0 = 3.0*y(2) - 2.0*y(3) call medof3(y0,y(1),y(2),ymed,change) y(1) = ymed c--- right end y0 = 3.0*y(n-1)-2.0*y(n-2) call medof3(y0,y(n),y(n-1),ymed,change) y(n) = ymed 999 return end subroutine split (y,n,change) c--- find 2-flats in y() and apply splitting algorithm integer n real y(n) logical change real w(6),y1 integer i1,i,nm2 c--- w() is a window 6 points wide which is slid along y() nm2 = n-2 do 10 i=1,4 w(i+2)= y(i) 10 continue c--- if y(1) = y(2) .ne. y(3), treat first two like a 2-flat with end-pt rule w(2) = y(3) i1 =1 20 if (w(3) .ne. w(4)) goto 40 if ( (w(3)-w(2)) * (w(5)-w(4)) .ge. 0.0) goto 40 c--- w(3) and w(4) form a two flat if (i1 .lt. 3) goto 30 c--- apply right end pt rule y1 = 3.0*w(2)-2.0*w(1) call medof3 (y1,w(3),w(2),y(i1),change) 30 if(i1 .ge. nm2) goto 40 c--- apply left end pt rule at i1+1 y1 = 3.0*w(5) - 2.0*w(6) call medof3 (y1,w(4),w(5),y(i1+1),change) c--- slide window 40 do 50 i=1,5 w(i) = w(i+1) 50 continue i1 = i1+1 if(i1 .ge. nm2) goto 60 w(6) = y(i1+3) goto 20 c--- apply rule to last 2 points if needed 60 w(6) = w(3) if(i1 .lt. n) goto 20 999 return end subroutine runmed (y,n,len,work,save,nw,err) c--- smooth y(0 by running medians of length len c note: use s3 for runnig medians of 3 instead of runmed integer n,len,nw,err real y(n),work(nw),save(nw) real median real temp,two integer savept,smopt,lenp1,i,j data two/2.0/ c--- work() is a local array in which data values are sorted c save acts as a window on the data if(len .le.nw) goto 5 err = 62 goto 999 5 do 10 i=1,len work(i) = y(i) save(i) = y(i) 10 continue savept = 1 smopt = int((float(len)+two)/two) lenp1 = len + 1 do 50 i=lenp1,n call qsort4 (len,work) if(err .ne. 0) goto 999 y(smopt) = median(work,len) temp = save(savept) do 20 j=1,len if(work(j) .eq. temp) goto 30 20 continue err = 63 goto 999 30 work(j) = y(i) save(savept) = y(i) savept = mod(savept,len)+1 smopt = smopt + 1 50 continue call qsort4 (len,work) if(err .ne. 0) goto 999 y(smopt) = median (work,len) 999 return end SUBROUTINE LOWESS (X,Y,N,F,NSTEPS,DELTA,YS,RW,RES) c--- LOWESS computes the smooth of a scatterplot of Y against X using c robust locally weighted regression. Fitted values, XS, are computed c at each of the values of the horizontal axis in X c--- see Cleveland, William S. (1979), J.A.S.A. 74, 829 c--- ARGUMENTS c X = input; abscissas of the points on the scatterplot; the values c of X must be ORDERED low to high on input c Y = input; ordinates of the points on the scatterplot c N = input; dimension of X, Y, YS, RW, and RES c F = input; specifies the amount of smoothing; F is the fraction c of the points used to compute each fitted value; as F increases c the smoothed values become smoother; choosing F in the c range 0.2 to 0.8 should serve most purposes; if you have no c idea what value to try, use F = 0.5 c NSTEPS = input; the number of iterations in the robust fit; if NSTEPS = 0, c the nonrobust fit is returned; setting NSTEPS = 2 should serve c most purposes c DELTA = input; non-negative parameter which may be used to save c computations; if N is less than 100, set DELTA = 0.0; if N c is greate than 100 you need to know some more about how DELTA c works c YS = output; fitted values; YS(I) is the fitted value at X(I); to c summarize the scatterplot, YS(I) should be plotted against c XS(I) c RW = output; robustness weights; RW(I) is the weight given to the c point (X(I),Y(I)); if NSTEPS = 0, RW is not used c RES = output; residuals; RES(I) = Y(I) - YS(I) c******************************************************************************* integer n integer nsteps real x(n),y(n),f,delta,ys(n),rw(n) real res(n) integer nright,min0,max0,i,j,ifix integer iter,last,m1,m2,ns,nleft real abs,cut,cmad,r,d1,d2 real c1,c9,alpha,denom,float logical ok if(n.ge.2) goto 1 ys(1) = y(1) return c--- at least 2, at most N points 1 ns = max0(min0(ifix(f*float(n)),n),2) iter = 1 goto 3 2 iter = iter+1 3 if (iter.gt.nsteps+1) goto 22 c--- robustness iterations nleft = 1 nright = ns c--- index of previous estimated point last = 0 c--- index of current point i = 1 4 if(nright.ge.n) goto 5 c--- move NLEFT, NRIGHT to right if radius decreases d1 = x(i)-x(nleft) c--- if D1 <= D2 with X(NRIGHT+1) == X(NRIGHT), LOWEST fixes d2 = x(nright+1)-x(i) if(d1.le.d2) goto 5 c--- radius will not decrease by move right nleft = nleft +1 nright = nright +1 goto 4 c--- fitted value at X(1) 5 CALL LOWEST (X,Y,N,X(I),YS(I),NLEFT,NRIGHT,RES,ITER + .GT.1,RW,OK) if (.not.ok) ys(i) = y(i) c--- all weights zero - copy over value (all RW = 0) if (last.ge.i-1) goto 9 denom = x(i)-x(last) c--- skipped points -- interpolate c--- non-zero - proof ? j = last+1 goto 7 6 j = j+1 7 if(j.ge.i) goto 8 alpha = (x(j)-x(last))/denom ys(j) = alpha*ys(i)+(1.0-alpha)*ys(last) goto 6 8 continue c--- last point actually estimated 9 last = i c--- X coord of close points cut = x(last)+delta i = last+1 goto 11 10 i = i+1 11 if (i.gt.n) goto 13 c--- find close points if (x(i).gt.cut) goto 13 c--- I one beyond last point within cut if (x(i).ne.x(last)) goto 12 ys(i) = ys(last) c--- exact match in X last = i 12 continue goto 10 c--- back 1 point so interpolation within DELTA, but always go forward 13 i = max0(last+1,i-1) 14 if (last.lt.n) goto 4 c--- residuals do 15 i=1,n res(i) = y(i) - ys(i) 15 continue if (iter.gt.nsteps) goto 22 c--- compute robustness weights except last time do 16 i=1,n rw(i) = abs(res(i)) 16 continue CALL QSORT4 (N,RW) m1 = n/2+1 m2 = n-m1+1 c--- median absolute residual cmad = 3.0*(rw(m1)+rw(m2)) c9 = 0.999*cmad c1 = 0.001*cmad do 21 i=1,n r = abs(res(i)) if (r.gt.c1) goto 17 rw(i) = 1. c--- near 0, avoid underflow goto 20 17 if (r.le.c9) goto 18 rw(i) = 0. c--- near 1, avoid underflow goto 19 18 rw(i) = (1.0-(r/cmad)**2)**2 19 continue 20 continue 21 continue goto 2 22 return end SUBROUTINE LOWEST (X,Y,N,XS,YS,NLEFT,NRIGHT,W,USERW,RW,OK) c--- LOWEST is a support routine for LOWESS and ordinarily will not be c called by the user. The fitted value, YS, is computed at the value c XS, of the horizontal axis. Robustness weights, RW, can be employed c in computing the fit. c--- see Cleveland, William S. (1979), J.A.S.A. 74, 829 c--- ARGUMENTS c X = input; abscissas of the points on the scatterplot; the values c of X must be ORDERED low to high on input c Y = input; ordinates of the points on the scatterplot c N = input; dimension of X, Y, YS, and RW c XS = input; value of the horizontal axis at which the smooth is c computed c YS = output; fitted value at XS c NLEFT = input; index of the first point which should be considered c in computing the fitted value c NRIGHT = input; index of the last point which should be considered c in computing the fitted value c W = output; W(I) is the weight for Y(I) used in the expression c for YS, which is the sum from I = NLEFT to NRIGHT of c W(I)*Y(I); W(I) is defined only a locations NLEFT to NRIGHT c USERW = input; logical variable; if USERW is >TRUE., a robust fit is c carried out using the weights in RW; if USERW is .FALSE., the c values in RW are not used c RW = input; robustness weights c OK = output; logical variable; if the weights for the smooth are c all 0.0, the fitted value YS is not computed and OK is set c to .FALSE.; if the fitted value is computed OK is set equal c to .TRUE. c******************************************************************************* integer n integer nleft,nright real x(n),y(n),xs,ys,w(n),rw(n) logical userw,ok integer nrt,j real abs,a,b,c,h,r real h1,sqrt,h9,amax1,range range = x(n) - x(1) h = amax1(xs-x(nleft),x(nright)-xs) h9 = 0.999*h h1 = 0.001*h c--- sum of weights a = 0.0 j = nleft goto 2 1 j = j+1 2 if(j.gt.n) goto 7 c--- compute weights (pick up all ties on right) w(j) = 0.0 r = abs(x(j)-xs) if (r.gt.h9) goto 5 if (r.le.h1) goto 3 w(j) = (1.0-(r/h)**3)**3 c--- small enough for non-zero weight goto 4 3 w(j) = 1.0 4 if (userw) w(j) = rw(j)*w(j) a = a+w(j) goto 6 5 if (x(j).gt.xs) goto 7 c--- get out at first zero at right 6 continue goto 1 c--- rightmost point (may be greater than NRIGHT because of ties) 7 nrt = j-1 if (a.gt.0.0) goto 8 ok = .false. goto 16 8 ok = .true. c--- weighted least squares do 9 j = nleft,nrt c--- make sum of W(J) = 1 w(j) = w(j)/a 9 continue if(h.le.0.0) goto 14 a = 0.0 c--- use linear fit do 10 j=nleft,nrt c--- weighted center of X values a = a+w(j)*x(j) 10 continue b = xs - a c = 0.0 do 11 j = nleft,nrt c = c+w(j)*(x(j)-a)**2 11 continue if(sqrt(c).le.0.001*range) goto 13 b = b/c c--- points are spread out enough to compute slope do 12 j = nleft,nrt w(j) = w(j)*(b*(x(j)-a)+1.0) 12 continue 13 continue 14 ys = 0.0 do 15 j = nleft,nrt ys = ys+w(j)*y(j) 15 continue 16 return end