SECTION 15.16 CONVOLUTIONS

15.16.1 CONTINUOUS RANDOM VARIABLES (CONVOLVE)

Given random variables x and y, with density functions f(x) and g(y) over the support (-infinity,+infinity), we compute the density function of the sum, x + y.

The MAIN program should include

        COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
        CALL DFLT

The convolution is computed by calling 
        SUBROUTINE CONVOLVE(Y,ANS,N,YL,YH,AX,BX,NPOINT,FUNC)

where 

        Y      = a DOUBLE PRECISION array dimensioned Y(N+1) which
                 will contain the abscissas at which the density is
                 calculated (output)
        ANS    = a DOUBLE PRECISION array dimensioned ANS(N+1) which
                 will contain the correponding densities (output)
        N      = the number of equal width subsegments into which
                 the range of the abscissas will be divided (input)
        YL,YH  = double precision scalars containing the lower and
                 upper bound of the abscissa range (input)
        AX,BX  = the lower and upper bound of abscissa values over
                 which numerical integration will be performed (input)
                 (DOUBLE PRECISION)
        NPOINT = the number of points used in iterated Gauss-Legendre
                 quadrature (input)
        FUNC   = the name of the function that evaluates the
                 integrand f(x)g(z-x), where z will be the
                 abscissa of the convolution (input). FUNC must
                 declared in an EXTERNAL statement in the MAIN program

The precise roles of YL, YH, AX and AB can be seen from noting that the convolution density at z is given by the integral

h(z)=IAXBXf(x)g(z-x)dx

Hence, for each point z, the numerical integration is from AX to BX; the convolution in turn is evaluated at N points along the z-scale between the values YL and YH. Thus, for example, the function HUNC in the following statements is the correct function for computing the convolution of two independent normals (not that one would ever want to use this method to compute that).


        function fdensty(x)
        implicit real*8 (a-h,o-z)
        pi=3.14159265
        sq2pi=dsqrt(2.d0*pi)
        xmu=0.
        sig=1.
        fdensty=1./(sq2pi*sig)*dexp(-0.5*(x-xmu)**2/sig**2)
        return
        end
        function gdensty(x)
        implicit real*8 (a-h,o-z)
        common/passy/z
        pi=3.14159265
        sq2pi=dsqrt(2.d0*pi)
        xmu=0.
        sig=10.
        gdensty=1./(sq2pi*sig)*dexp(-0.5*(z-x-xmu)**2/sig**2)
        return
        end
        function hunc(x)
        implicit real*8 (a-h,o-z)
        hunc=fdensty(x)*gdensty(x)
        return
        end
Notice that either HUNC or other routines it calls (in this case GDENSTY(X) must contain the statement COMMON/PASSY/Z, so that the subroutine CONVOLVE is able to pass different Z values to the integrand evaluator.

Care has to be exercized so that the entire relevant range of the convolution is included and that the integrations that are performed (one for each value of Y) cover the relevant range. For example, for the problem above, it is reasonable to set AX=-60., BX=60., NPOINT=160, N=500, YL=-30., YH=30. If, after the computation, the area under the computed density is not sufficiently close to 1.0, some of these parameters have to be reset.

Return to

|Beginning|SHELL|

15.16.2 DISCRETE RANDOM VARIABLES (DISCCONV)

Given discrete probability distributions f(x), x=0,1,2,...,n and g(y), y=0,1,2,...,m, DISCCONV computes the discrete probability distribution of z = x + y, z=0,1,2,...,m+n.

The MAIN program should include

        COMMON/BPRINT/IPT,NFILE,NDIG,NPUNCH,JPT,MFILE
        CALL DFLT

The convolution is computed by 
      CALL DISCCONV(F,G,N,M,ANS,SUM,XMEAN,STDVAR,FUNCF,FUNCG)
where 
        F      = a DOUBLE PRECISION array dimensioned F(0:N) which
                 will contain the probabilities of the x variable
                 (output)
        G      = a DOUBLE PRECISION array dimensioned G(0:M) which
                 will contain the probabilities of the y variable
                 (output)
        N      = number of outcomes for x minus 1 (input)
        M      = number of outcomes for y minus 1 (input)
        ANS    = a DOUBLE PRECISION array dimensioned ANS(0:N+M)
                 (output)
        SUM    = DOUBLE PRECISION scalar, will contain the sum of
                 the probabilities in ANS (output)
        XMEAN  = DOUBLE PRECISION scalar, will contain mean value
                 of x+y (output)
        STDVAR = DOUBLE PRECISION scalar, will contain standard
                 deviation of x+y (output)
        FUNCF  = name of DOUBLE PRECISION function that evaluates
                 the probabilities of x (input, has to be declared
                 in an EXTERNAL statement and user must provide
                 the program for it)
        FUNCG  = name of DOUBLE PRECISION function that evaluates
                 the probabilities of y (input, has to be declared
                 in an EXTERNAL statement and user must provide
                 the program for it)
If one or both of the random variables can attain infintely many values, an approximation to the convolution may still be obtained by setting N and or M very large. The quality of the approximation will depend on how close SUM is to 1.0.
Return to

|Beginning|SHELL|