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 DFLTThe convolution is computed by calling
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
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 endNotice 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.
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 DFLTThe convolution is computed by
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.