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 ) C 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 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