C     PROGRAM 4.2  KLINFO
C
C  ...  Driver program of the subroutine KLINFO  ...
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  PARAMG(3), PARAMF(3)
      EXTERNAL  GAUSS
      EXTERNAL  CAUCHY
      DATA  PARAMG /0.0D0, 1.0D0, 0.0D0/
      DATA  PARAMF /0.1D0, 1.5D0, 0.0D0/
C
      XMIN = -8.0D0
      XMAX =  8.0D0
      WRITE(6,600)  XMIN, XMAX
      WRITE(6,610)
      DO 10 II=1,4
      NINT = (XMAX-XMIN+1.0D-5)*2**(II-1)
   10 CALL KLINFO(GAUSS,GAUSS ,PARAMG,PARAMF,XMIN,XMAX,NINT,FKLI,GINT)
      STOP
  600 FORMAT( 1H ,'XMIN =',F5.1,3X,'XMAX =',F5.1 )
  610 FORMAT( 1H ,'  XMIN  XMAX   NINT',3X,'DX',9X,'FKLI',11X,'GINT' )
      E N D
      SUBROUTINE  KLINFO( DISTG,DISTF,PARAMG,PARAMF,XMIN,XMAX,NINT,
     *                    FKLI,GINT )
C
C  ...  This subroutine computes Kullback-Leibler information  ...
C
C     Inputs:
C        DISTG:    function name for the true density
C        DISTF:    function name for the model density
C        PARAMG:   parameter vector of the true density
C        PARAMF:   parameter vector of the model density
C        XMIN:     lower limit of integration
C        XMAX:     upper limit of integration
C        NINT:     number of function evaluation
C     Outputs:
C        FKLI:     Kullback-Leibler information number, I(g;f)
C        GINT:     integration of g(y) over [XMIN,XMAX]
C
      IMPLICIT REAL*8(A-H,O-Z)
      DIMENSION  PARAMG(*), PARAMF(*)
C
      DX = (XMAX-XMIN)/DFLOAT(NINT)
      FKLI = 0.0D0
      GINT = 0.0D0
C
      DO 10  I=0,NINT
      XX = XMIN + DX*I
      GX = DISTG( XX,PARAMG )
      FX = DISTF( XX,PARAMF )
      IF( I.GE.1 .AND. I.LT.NINT )  THEN
         FKLI = FKLI + DLOG( GX/FX )*GX
         GINT = GINT + GX
      ELSE
         FKLI = FKLI + DLOG( GX/FX )*GX/2.0D0
         GINT = GINT + GX/2.0D0
      END IF
   10 CONTINUE
C
      FKLI = FKLI*DX
      GINT = GINT*DX
      WRITE(6,600) XMIN, XMAX, NINT, DX, FKLI, GINT
C
      RETURN
  600 FORMAT( 1H ,2F6.2,I5,F9.4,D17.8,F12.8 )
      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