C PROGRAM 4.1 DENSTY C C ... This program draws probability density function ... C C The following inputs are required in the main program: C MODEL: function number (1: GAUSS, 2: CAUCHY, etc.) C XMIN: lower bound of the interval C XMAX: upper bound of the interval C PARAM: parameter vector C @TEST.PN41: DEC.26,1990, 8/6/91 C PARAMETER( K=201 ) IMPLICIT REAL*8(A-H,O-Z) CHARACTER*24 TITLE1(0:7) CHARACTER*72 TITLE DIMENSION F(K), PARAM(3), NP(0:7) COMMON /CMDATA/ TITLE EXTERNAL USERF EXTERNAL GAUSS EXTERNAL CAUCHY EXTERNAL PEARSN EXTERNAL EXPNTL EXTERNAL CHISQR EXTERNAL DBLEXP EXTERNAL UNIFRM C DATA TITLE1/'USER SUPPLIED FUNCTION ','NORMAL DISTRIBUTION ' * ,'CAUCHY DISTRIBUTION ','PEARSON DISTRIBUTION ' * ,'EXPONENTIAL DISTRIBUTION','CHI-SQUARE DISTRIBUTION ' * ,'DOUBLE EXPONENTIAL DIST.','UNIFORM DISTRIBUTION '/ DATA NP/0,2,2,3,1,1,0,2/ C WRITE(6,*) 'INPUT MODEL NUMBER ?' READ(5,*) MODEL WRITE(6,*) 'INPUT XMIN AND XMAX ?' READ(5,*) XMIN, XMAX IF( MODEL.EQ.0 ) THEN WRITE(6,*) 'INPUT THE NUMBER OF PARAMETERS' READ(5,*) NP(0) END IF WRITE(6,*) 'INPUT PARAM(I),I=1,',NP(MODEL),' ?' READ(5,*) (PARAM(I),I=1,NP(MODEL)) TITLE = TITLE1(MODEL) DO 10 I=1,6 10 TITLE = TITLE//' ' C IF(MODEL.EQ.1) CALL DENSTY( GAUSS ,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.2) CALL DENSTY( CAUCHY,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.3) CALL DENSTY( PEARSN,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.4) CALL DENSTY( EXPNTL,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.5) CALL DENSTY( CHISQR,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.6) CALL DENSTY( DBLEXP,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.7) CALL DENSTY( UNIFRM,F,K,PARAM,XMIN,XMAX ) IF(MODEL.EQ.0) CALL DENSTY( USERF ,F,K,PARAM,XMIN,XMAX ) CALL PTDENS( F,K,XMIN,XMAX,PARAM,NP(MODEL) ) CALL PRDENS( F,K ) STOP 600 FORMAT( 1H ,10F8.4 ) E N D SUBROUTINE PRDENS( F,K ) C C ... This subroutine prints out density ... C C Inputs: C F(I): values of density function C K: data length C IMPLICIT REAL*8(A-H,O-Z) DIMENSION F(K) C WRITE(6,600) WRITE(6,610) (F(I),I=1,K) C RETURN 600 FORMAT( 1H ,'PROGRAM 4.1: ' ) 610 FORMAT( 1H ,10F8.4 ) E N D SUBROUTINE PTDENS( F,K,XMIN,XMAX,PARAM,NP ) C C ... This subroutine draws density function ... C C Inputs: C F(I): values of density function C K: data length C XMIN: the lower edge of the interval C XMAX: the upper edge of the interval C PARAM: parameter vector C NP: number of parameters C IMPLICIT REAL*8(A-H,O-Z) CHARACTER VNAME*8 DIMENSION F(K), PARAM(3), VNAME(5) C WX = 12.0 WY = 12.0 VNAME(1) = 'PARAM(1)' VNAME(2) = 'PARAM(2)' VNAME(3) = 'PARAM(3)' C CALL PLOTS C call plots( 1,0,0,1,0 ) C call form( 1 ) C call factor( 10.0 ) CALL HEADER( 'PROGRAM 4.1: DENSITY FUNCTION',30,NP,VNAME,PARAM) CALL DELX( XMIN,XMAX,DX ) CALL MAXMIN( F,K,YMIN,YMAX,DY ) CALL DRAWY2( 'DENSITY FUNCTION',16,3.0D0,15.0-WY,F,K,WX,WY, * 0,0,XMIN,XMAX,DX,YMIN,YMAX,DY ) CALL PLOTE C call plot( 0.0,0.0,999 ) C RETURN E N D SUBROUTINE DENSTY( DIST,P,K,PARAM,XMIN,XMAX ) C C ... This subroutine evaluates values of density ... C DIST(X), X=XMIN,...,XMAX C Inputs: C DIST: name of function C PARAM: parameters of the density C XMIN: minimum of X C XMAX: maximum of X C K: number of location, I-th location is XMIN + I*DX C where DX = (I-1)*(XMAX-XMIN)/(K-1) C OUTPUT: C P(I): density of DIST at I-th location C IMPLICIT REAL*8( A-H,O-Z ) DIMENSION P(K), PARAM(*) EXTERNAL DIST C DX = (XMAX-XMIN)/(K-1) DO 10 I=1,K X = XMIN + DX*(I-1) 10 P(I) = DIST( X,PARAM ) RETURN E N D DOUBLE PRECISION FUNCTION GAUSS( X,PARAM ) C C ... Gaussian (normal) distribution ... C C Inputs: C X: C PARAM(1): mean C PARAM(2): variance C Output: C GAUSS: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) DATA C1 /2.506628275D0/ C GAUSS = DEXP( -(X-PARAM(1))**2/(2*PARAM(2)) )/(C1*DSQRT(PARAM(2))) RETURN E N D DOUBLE PRECISION FUNCTION CAUCHY( X,PARAM ) C C ... Cauchy distribution ... C C Inputs: C X: C PARAM(1): location parameter, mu C PARAM(2): dispersion parameter, tau square C Output: C CAUCHY: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) DATA PI /3.1415926535D0/ C CAUCHY = DSQRT( PARAM(2) )/(PARAM(2) + (X-PARAM(1))**2)/PI RETURN C E N D DOUBLE PRECISION FUNCTION PEARSN( X,PARAM ) C C ... Pearson family of distributions ... C C Inputs: C X: C PARAM(1): location parameter, mu C PARAM(2): dispersion parameter, tau square C PARAM(3): shape parameter C Output: C PEARSN: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(3) DATA PI/3.1415926535D0/ C PEARSN = DGAMMA(PARAM(3))/DGAMMA(PARAM(3)-0.5D0) * /DSQRT(PI)*PARAM(2)**(PARAM(3)-0.5D0) * /((X-PARAM(1))**2 + PARAM(2))**PARAM(3) RETURN C END DOUBLE PRECISION FUNCTION DBLEXP( X,PARAM ) C C ... double exponential distribution f(x) = exp(x - exp(x)) ... C C Inputs: C X: C Output: C DBLEXP: density at X C IMPLICIT REAL*8(A-H,O-Z) C DBLEXP = DEXP( X-DEXP(X) ) RETURN C E N D DOUBLE PRECISION FUNCTION EXPNTL( X,PARAM ) C C ... Exponential distribution ... C C Inputs: C X: C PARAM(1): lambda C Output: C EXPNTL: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(1) C IF( X.GE.0.0D0 ) EXPNTL = PARAM(1)*DEXP( -PARAM(1)*X ) IF( X.LT.0.0D0 ) EXPNTL = 0.0D0 RETURN E N D DOUBLE PRECISION FUNCTION CHISQR( X,PARAM ) C C ... Chi-square distributions ... C C Inputs: C X: C PARAM(1): degree of freedoms, k C Output: C CHISQR: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(*) C IF( X.GT.0.0D0 ) CHISQR = DEXP( -X/2 )*(X/2)**(PARAM(1)/2-1.D0) * /(2*DGAMMA(PARAM(1)/2)) IF( X.LE.0.0D0 ) CHISQR = 0.0D0 RETURN E N D DOUBLE PRECISION FUNCTION UNIFRM( X,PARAM ) C C ... Uniform distribution on [a,b] ... C C Inputs: C X: C PARAM(1): a C PARAM(2): b C Output: C UNIFRM: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) C IF( X.GT.PARAM(1) .AND. X.LE.PARAM(2) ) THEN UNIFRM = 1.0D0/(PARAM(2)-PARAM(1)) ELSE UNIFRM = 0.0D0 END IF RETURN E N D DOUBLE PRECISION FUNCTION USERF( X,PARAM ) C C ... User supplied density function ... C (The following is an example of two-sided exponential dist.) C C Inputs: C X: C PARAM(1): lambda C Output: C USERF: density at X C IMPLICIT REAL*8(A-H,O-Z) DIMENSION PARAM(2) C IF( X.GE.0.0D0 ) THEN USERF = PARAM(1)*DEXP( -PARAM(1)*X )/2 ELSE USERF = PARAM(1)*DEXP( PARAM(1)*X )/2 END IF RETURN E N D DOUBLE PRECISION FUNCTION DGAMMA( X ) C C ... Gamma function ... C IMPLICIT REAL*8(A-H,O-Z) DIMENSION A(0:10) DATA A /0.999999999999269D0, 0.42278433696202D0, * 0.41184025179616D0, 0.08157821878492D0, 0.0742379070629D0, * -0.00021090746731D0, 0.01097369584174D0,-0.00246674798054D0, * 0.00153976810472D0,-0.00034423420456D0, 0.00006771057117D0/ C DGAM = 1.0D0 Y = X IF( X.GT.3.0D0 ) THEN 10 Y = Y-1 DGAM = DGAM*Y IF( Y.GT.3.0D0 ) GO TO 10 END IF IF( X.LE.2.0D0 ) THEN 20 DGAM = DGAM/Y Y = Y+1 IF( Y.LE.2.0D0 ) GO TO 20 END IF C Z = 1.0D0 SUM = 0.0D0 DO 30 I=0,10 SUM = SUM + A(I)*Z 30 Z = Z*(Y-2) DGAMMA = DGAM*SUM RETURN E N D SUBROUTINE AXISXY(X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC,IBOX,IX,IY) C C ... This subroutine draws X and Y axes. C C Inputs: C X,Y: location of the left bottom of the figure C WX,WY: width and height of the figure C X0,X1: lower and upper bounds of the X axis C Y0,Y1: lower and upper bounds of the Y axis C DX,DY: step width in X and Y axes C DWC: size of characters C IBOX: = 0 draw X and Y axes only C = 1 draw window doundary C IX,IY: number of additional click in each step C Modified: 8/31/90, 11/21/90 C REAL*8 X,Y,WX,WY,X0,X1,Y0,Y1,DX,DY,DWC C FX = X FY = Y FWX = WX FWY = WY FWC = DWC CALL NEWPEN( 2 ) CALL PLOT( FX,FY,-3 ) CALL PLOT( 0.0,FWY,3 ) CALL PLOT( 0.0,0.0,2 ) CALL PLOT( FWX,0.0,2 ) IF ( IBOX.EQ.1 ) THEN CALL PLOT( FWX,FWY,2 ) CALL PLOT( 0.0,FWY,2 ) END IF C C ... draw X axis ... C WC = DWC NX = (X1-X0)/DX + 0.001 NY = (Y1-Y0)/DY + 0.001 MX = -DLOG10( DX ) + 0.001 IF(DX.LT.1.0D0) MX = -DLOG10( DX ) + 0.999 IF( MX.LE.0 ) MX = -1 D = WX*DX/(X1-X0) DO 10 I=1,NX+1 XN = X0 + DX*(I-1) XX = D*(I-1) CALL PLOT( XX,-FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF LX = 0 IF(XN.GT.0.0D0) LX = ALOG10( XN ) + 0.001 IF( MX.GT.0 ) XXX = WC*(LX+MX)*0.5 + WC*4/7.0 IF( MX.LT.0 ) XXX = WC*LX*0.5 + WC*4/7.0 CALL NUMBER( XX-XXX,-(1.5*FWC+0.1),FWC,XN,0.0,MX ) 10 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 20 I=1,NX+1 DO 20 J=1,IX XX = D*(I-1) + (D*J)/IX IF( XX.LE.WX ) THEN CALL PLOT( XX,FWC/2,3 ) CALL PLOT( XX,0.0,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( XX,FWY-FWC/2,3 ) CALL PLOT( XX,FWY,2 ) END IF END IF 20 CONTINUE CALL NEWPEN( 2 ) END IF C C ... draw Y axis ... C D = WY*DY/(Y1-Y0) MY = -DLOG10( DY ) + 0.99999 IF( MY.LE.0 ) MY = -1 DO 30 I=1,NY+1 YN = Y0 + DY*(I-1) YY = D*(I-1) CALL PLOT( -FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF LY = 1 IF(YN.LT. 0.0) LY = 2 IF(YN.GE. 1.0) LY = ALOG10( YN )+1.00001 IF(YN.LE.-1.0) LY = ALOG10(-YN )+2.00001 IF( MY.GT.0 ) YYY = WC*(LY+MY+1) + 0.20 IF( MY.LT.0 ) YYY = WC*LY + 0.20 CALL NUMBER( -YYY,YY-FWC/2,FWC,YN,0.0,MY ) 30 CONTINUE IF( IX.GT.1 ) THEN CALL NEWPEN(1) DO 40 I=1,NY+1 DO 40 J=1,IY YY = D*(I-1) + (D*J)/IY IF( YY.LE.WY ) THEN CALL PLOT( FWC/2,YY,3 ) CALL PLOT( 0.0,YY,2 ) IF( IBOX.EQ.1 ) THEN CALL PLOT( FWX-FWC/2,YY,3 ) CALL PLOT( FWX,YY,2 ) END IF END IF 40 CONTINUE CALL NEWPEN( 2 ) END IF C RETURN E N D SUBROUTINE HEADER( PNAME,NP,M,VNAME,VALUE ) C C ... This subroutine draw title of data, program name, names C and values of parameters ... C C Inputs: C PNAME: program name C NP: number of characters of program name C M: number of parameters C VNAME(I): name of the I-th parameter C VALUE(I): value of the I-th parameter C Input via common: C TITLE: title of the data set C CHARACTER PNAME*80, TITLE*72 CHARACTER*8 VNAME(M) REAL*8 VALUE DIMENSION VALUE(M), DAY(2), TIME(3) COMMON /CMDATA/ TITLE C CALL NEWPEN( 2 ) CALL DATE( DAY ) CALL CLOCK( TIME,1 ) C CALL SYMBOL( 2.0,18.7,0.25,TITLE,0.0,72 ) CALL SYMBOL( 2.0,18.3,0.20,PNAME,0.0,NP ) CALL NEWPEN( 1 ) CALL SYMBOL( 2.0,17.9,0.20,DAY,0.0,8 ) CALL SYMBOL( 4.0,17.9,0.20,TIME,0.0,12 ) YY = 17.9 DO 10 I=1,M IF( MOD(I,4).EQ.1 ) THEN YY = YY - 0.40 X1 = 2.0 X2 = 3.5 ELSE X1 = X1 + 5.0 X2 = X2 + 5.0 END IF CALL SYMBOL( X1,YY,0.20,VNAME(I),0.0,8 ) 10 CALL NUMBER( X2,YY,0.20,SNGL(VALUE(I)),0.0,6 ) C RETURN E N D SUBROUTINE DELX( XMIN,XMAX,DX ) C C ... This subroutine determines step width in drawing axis ... C C Inputs: C XMIN: the minimum value in the axis C XMAX: the maximum value in the axis C Output: C DX: step size C IMPLICIT REAL*8(A-H,O-Z) DDX = XMAX - XMIN DX = INT( DLOG10( DDX ) ) IF( DDX.LT.1.0D0 ) DX = DX - 1 DX = 10.0D0**DX IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 IF( DDX/DX.LT.3.0D0 ) DX = DX/2.0 RETURN E N D SUBROUTINE DRAWY2( FNAME,NC,XO,YO,Y,N,WX,WY,IPOS,ISW,XMIN, * XMAX,DX,YMIN,YMAX,DY ) C C ... This subroutine draws time series data ... C C Inputs: C FNAME: title of the data set C NC: number of charcters of the name C XO,YO: origin of the figure C Y(I): time series C N: data length C WX,WY: width and height of the figure C IPOS: starting position of time series C ISW: = 1 connect data by line C = 2 connect data by dashed line C XMIN,XMAX: minimum and maximum of the X axis C DX: step width of the X axis C YMIN,YMAX: minimum and maximum of the Y axis C DY: step width of the Y axis C IMPLICIT REAL*8(A-H,O-Z) CHARACTER FNAME*80 DIMENSION Y(N) C CALL AXISXY(XO,YO,WX,WY,XMIN,XMAX,YMIN,YMAX,DX,DY,0.25D0,1,10,1) IF( ISW.LE.1 ) CALL PLOTY ( Y,N,YMIN,YMAX,WX,WY,IPOS,1 ) IF( ISW.EQ.1 ) CALL PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,0.0D0 ) IF( ISW.EQ.2 ) CALL PLOTY2( Y,N,YMIN,YMAX,WX,WY,IPOS,1,YMIN ) CALL NEWPEN( 1 ) CALL SYMBOL( 0.0,SNGL(WY)+0.2,0.3,FNAME,0.0,NC ) RETURN E N D SUBROUTINE MAXMIN( X,N,XMIN0,XMAX0,DXL ) C C ... This subroutine determines the minimum, the maximum and C the step width ... C C Inputs: C X(I): data C N: data length C Outputs: C XMIN0: the minimum bound for the figure C XMAX0: the maximum bound for the figure C DXL: step width C IMPLICIT REAL*8(A-H,O-Z) DIMENSION X(N) C XMIN = 1.0D30 XMAX =-1.0D30 C DO 10 I=1,N IF( X(I) .LT. XMIN ) XMIN = X(I) 10 IF( X(I) .GT. XMAX ) XMAX = X(I) DX = XMAX-XMIN IF( DLOG10(DX) .GE. 0.0D0 ) DXL = INT( DLOG10(DX) ) IF( DLOG10(DX) .LT. 0.0D0 ) DXL = INT( DLOG10(DX) )-1.0 DXL = 10.0**DXL IF( DX/DXL.GT.6.0D0 ) DXL = DXL*2.0 DIF = INT( DX/DXL ) XMIN0 = INT( XMIN/DXL )*DXL XMAX0 = XMIN0 + DIF*DXL IF( XMIN0 .GT. XMIN ) XMIN0 = XMIN0 - DXL 30 IF( XMAX0 .GE. XMAX ) GO TO 40 XMAX0 = XMAX0 + DXL GO TO 30 40 CONTINUE C RETURN E N D SUBROUTINE PLOTY( Y,N,YMIN,YMAX,WX,WY,IPOS,ISW ) C C ... This subroutine draws time vs data plot of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 connect data by straight line C = 2 connect data by dashed line C = 3 express data by character (character code=ISW) C Modified: 8/31/90 C REAL*8 Y(N), YMIN, YMAX, WX, WY C DX = WX/(N-1+IPOS) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IPOS) YY = (Y(I) - YMIN)*DY IF( I.EQ.1 .AND. ISW.LT.3 ) THEN CALL PLOT( XX,YY,3 ) ELSE IF(ISW.EQ.1) CALL PLOT( XX,YY,2 ) IF(ISW.EQ.2) CALL DASHP( XX,YY,0.2 ) IF(ISW.GE.3) CALL SYMBOL( XX,YY,0.2,ISW,0.0,-1 ) END IF 100 CONTINUE C RETURN E N D SUBROUTINE PLOTY2( Y,N,YMIN,YMAX,WX,WY,IOFF,ISW,YLEVEL ) C C ... This subroutine draws box car graph of the data ... C C Inputs: C Y(I): data C N: data length C YMIN: minimum bound for the Y axis C YMAX: maximum bound for the Y axis C WX,WY: width and height of the figure C IPOS: starting position of the data C ISW: = 1 draw box with thick line C = 2 draw box with thin line C YLEVEL: ground level for the box car C REAL*8 Y(N), YMIN, YMAX, WX, WY, YLEVEL C DX = WX/(N-1+IOFF) DY = WY/(YMAX - YMIN) C CALL NEWPEN( 3 ) IF( N.GT.75 ) CALL NEWPEN( 2 ) IF( N.GT.150) CALL NEWPEN( 1 ) IF( ISW.EQ.2 ) CALL NEWPEN( 1 ) DO 100 I=1,N XX = DX*(I-1+IOFF) YY = (Y(I) - YMIN)*DY Y0 = (YLEVEL-YMIN)*DY CALL PLOT( XX,Y0,3 ) CALL PLOT( XX,YY,2 ) 100 CONTINUE CALL NEWPEN( 2 ) C RETURN E N D