13 Processing the variance array

Errors associated with data can be accommodated in an NDF by means of a Variance array. Such an array has the same shape and size as the associated data array, and contains an estimate of the variance for each element in that data array. Note that it is the variance for each data element which is stored rather than the standard deviation, the former being simply the square of the latter. Variance is the chosen method of storage for errors as most error processing is considerably simpler when variance rather than standard deviation is used. For example if two numbers are added: z = x + y

the standard deviation in the result, σz, is related to the standard deviations in the input numbers σx, σy, as follows: σz = σx 2 + σy 2

whereas if the variance is considered, the variance in the result Vz is simply the sum of the variances of the input numbers, i.e. Vz = Vx + Vy

The latter is significantly quicker to compute.

Application programs should consider whether an input variance array remains valid after a data array has been processed. For example, the variance array remains valid when a constant with no associated error is added to each element of the main data array. In such a case the variance array should be propagated unchanged to the output NDF.

The program discussed in the previous section took the square root of each element in the main data array of the input NDF. Note that variance was not among the components included in the call to NDF_PROP, so any variance array in the input NDF would not be propagated to the output NDF. Obviously in this case, the input variance array is not appropriate to the output data, as taking the square root of a data set changes the variance.

However, in the case where Gaussian statistics apply and a data element is raised to a power, i.e. y = xn

the variance in y, Vy is related to the variance in x, Vx, as follows: Vy y2 = n2 Vx x2

In the case of taking a square root (n = 1/2) the variance is given by: Vy = Vx 4x

ADAM_EXAMPLES:SQROOTV.FOR processes the variance array according to this formula. In the interests of simplicity, the program does not try to cope with bad data values, but as described in the previous section it should ensure that it does not attempt the processing of data arrays which it cannot handle. (However in the form shown here, it will crash if the input data contain negative numbers! A more robust program is presented in the next section.)

The check for the presence of bad pixels in the input data is as follows:

  *  Abort if main data array contains bad pixels.
        CALL NDF_BAD (NDF1, ’Data’, .TRUE., BAD, STATUS)
        IF (BAD) THEN
           STATUS = SAI__ERROR
           CALL ERR_REP (’SQROOTV_BADPIX’,
       :                ’Sorry, cannot cope with bad pixels’, STATUS)
           GOTO 999
        ENDIF

The remainder of the program is summarised in the following steps. The input NDF is checked for the existence of a variance array; if such an array exists, then both the main data and variance arrays will be processed by the program, otherwise only the main data array will be processed.

NDF_MAP can be given a list of components rather than just one. In the example below the list will comprise either ’Data’ or ’Data,Variance’ as appropriate.11 When a list is supplied to NDF_MAP, the pointer argument must be an array of at least the same size as the number of components in the list. A pointer to each mapped component will be returned via the array, in the order corresponding to that in the component list.

The output NDF is created by the propagation of the input NDF. The variance is not among the components propagated as it is more efficient in this case to create a new structure in the output NDF rather than copy the existing variance array and overwrite it. The NDF_MAP call which maps the variance for ’WRITE’ in the output NDF creates a variance structure of the same size and type as the data array. If a different type is desired, the data and variance can be mapped separately with the chosen types, (of course the size of the variance array must match the data array).

The mapped arrays are passed to a subroutine which takes the square root of each element in the main data array and calculates the variance appropriate to the output NDF if a variance structure was found in the input NDF.

  *   Check whether variance array exists
        CALL NDF_STATE (NDF1, ’Variance’, VARNCE, STATUS)
  
  *   Work out component list for NDF_MAP. If a variance array exists,
  *   both data and variance are mapped - otherwise only the data.
        COMP = ’ ’
        IF (VARNCE) THEN
           COMP = ’Data, Variance’
        ELSE
            COMP = ’Data’
        ENDIF
  
  *  Create a new output NDF based on the input NDF.
        CALL NDF_PROP (NDF1, ’Axis’, ’OUTPUT’, NDF2, STATUS )
  
  *  Map the input and output data and variance arrays.
        CALL NDF_MAP (NDF1, COMP, ’_REAL’, ’READ’, PNTR1, NELM, STATUS)
        CALL NDF_MAP (NDF2, COMP, ’_REAL’, ’WRITE’, PNTR2, NELM, STATUS)
  
  *   Take the square root of each element in the main data array.
  *   The output variance is also generated if appropriate.
        CALL SQRTVR(NELM, VARNCE, %VAL(PNTR1(1)), %VAL(PNTR1(2)),
       :            %VAL(PNTR2(1)), %VAL(PNTR2(2)), STATUS)

And the code from the subroutine:

        SUBROUTINE SQRTVAR (NELM, VARNCE, IN, VARIN, OUT, VAROUT, STATUS)
  *     ...
  *   Take the square root of each input element and find the variance.
        DO I = 1, NELM
           OUT(I) = SQRT(IN(I))
           IF (VARNCE) VAROUT(I) = VARIN(I)*0.25/IN(I)
        ENDDO
        END

11One advantage of mapping the data and variance arrays with a single call is that it is more efficient for NDF_MAP to access any quality array only once and insert bad pixels into the mapped data and variance arrays simultaneously.