NAG Library Routine Document

g03dbf (discrim_mahal)

1
Purpose

g03dbf computes Mahalanobis squared distances for group or pooled variance-covariance matrices. It is intended for use after g03daf.

2
Specification

Fortran Interface
Subroutine g03dbf ( equal, mode, nvar, ng, gmn, ldgmn, gc, nobs, m, isx, x, ldx, d, ldd, wk, ifail)
Integer, Intent (In):: nvar, ng, ldgmn, nobs, m, isx(*), ldx, ldd
Integer, Intent (Inout):: ifail
Real (Kind=nag_wp), Intent (In):: gmn(ldgmn,nvar), gc((ng+1)*nvar*(nvar+1)/2), x(ldx,*)
Real (Kind=nag_wp), Intent (Inout):: d(ldd,ng)
Real (Kind=nag_wp), Intent (Out):: wk(2*nvar)
Character (1), Intent (In):: equal, mode
C Header Interface
#include <nagmk26.h>
void  g03dbf_ (const char *equal, const char *mode, const Integer *nvar, const Integer *ng, const double gmn[], const Integer *ldgmn, const double gc[], const Integer *nobs, const Integer *m, const Integer isx[], const double x[], const Integer *ldx, double d[], const Integer *ldd, double wk[], Integer *ifail, const Charlen length_equal, const Charlen length_mode)

3
Description

Consider p variables observed on ng populations or groups. Let x-j be the sample mean and Sj the within-group variance-covariance matrix for the jth group and let xk be the kth sample point in a dataset. A measure of the distance of the point from the jth population or group is given by the Mahalanobis distance, Dkj:
Dkj2=xk-x-jTSj-1xk-x-j.  
If the pooled estimated of the variance-covariance matrix S is used rather than the within-group variance-covariance matrices, then the distance is:
Dkj2=xk-x-jTS-1xk-x-j.  
Instead of using the variance-covariance matrices S and Sj, g03dbf uses the upper triangular matrices R and Rj supplied by g03daf such that S=RTR and Sj=RjTRj. Dkj2 can then be calculated as zTz where Rjz=xk-x-j or Rz=xk-x-j as appropriate.
A particular case is when the distance between the group or population means is to be estimated. The Mahalanobis squared distance between the ith and jth groups is:
Dij2=x-i-x-jTSj-1x-i-x-j  
or
Dij2=x-i-x-jTS-1x-i-x-j.  
Note:  Djj2=0 and that in the case when the pooled variance-covariance matrix is used Dij2=Dji2 so in this case only the lower triangular values of Dij2, i>j, are computed.

4
References

Aitchison J and Dunsmore I R (1975) Statistical Prediction Analysis Cambridge
Kendall M G and Stuart A (1976) The Advanced Theory of Statistics (Volume 3) (3rd Edition) Griffin
Krzanowski W J (1990) Principles of Multivariate Analysis Oxford University Press

5
Arguments

1:     equal – Character(1)Input
On entry: indicates whether or not the within-group variance-covariance matrices are assumed to be equal and the pooled variance-covariance matrix used.
equal='E'
The within-group variance-covariance matrices are assumed equal and the matrix R stored in the first pp+1/2 elements of gc is used.
equal='U'
The within-group variance-covariance matrices are assumed to be unequal and the matrices Rj, for j=1,2,,ng, stored in the remainder of gc are used.
Constraint: equal='E' or 'U'.
2:     mode – Character(1)Input
On entry: indicates whether distances from sample points are to be calculated or distances between the group means.
mode='S'
The distances between the sample points given in x and the group means are calculated.
mode='M'
The distances between the group means will be calculated.
Constraint: mode='M' or 'S'.
3:     nvar – IntegerInput
On entry: p, the number of variables in the variance-covariance matrices as specified to g03daf.
Constraint: nvar1.
4:     ng – IntegerInput
On entry: the number of groups, ng.
Constraint: ng2.
5:     gmnldgmnnvar – Real (Kind=nag_wp) arrayInput
On entry: the jth row of gmn contains the means of the p selected variables for the jth group, for j=1,2,,ng. These are returned by g03daf.
6:     ldgmn – IntegerInput
On entry: the first dimension of the array gmn as declared in the (sub)program from which g03dbf is called.
Constraint: ldgmnng.
7:     gcng+1×nvar×nvar+1/2 – Real (Kind=nag_wp) arrayInput
On entry: the first pp+1/2 elements of gc should contain the upper triangular matrix R and the next ng blocks of pp+1/2 elements should contain the upper triangular matrices Rj. All matrices must be stored packed by column. These matrices are returned by g03daf. If equal='E' only the first pp+1/2 elements are referenced, if equal='U' only the elements pp+1/2+1 to ng+1pp+1/2 are referenced.
Constraints:
  • if equal='E', R0.0;
  • if equal='U', the diagonal elements of the Rj0.0, for j=1,2,,ng.
8:     nobs – IntegerInput
On entry: if mode='S', the number of sample points in x for which distances are to be calculated.
If mode='M', nobs is not referenced.
Constraint: if nobs1, mode='S'.
9:     m – IntegerInput
On entry: if mode='S', the number of variables in the data array x.
If mode='M', m is not referenced.
Constraint: if mnvar, mode='S'.
10:   isx* – Integer arrayInput
Note: the dimension of the array isx must be at least max1,m.
On entry: if mode='S', isxl indicates if the lth variable in x is to be included in the distance calculations. If isxl>0 the lth variable is included, for l=1,2,,m; otherwise the lth variable is not referenced.
If mode='M', isx is not referenced.
Constraint: if mode='S', isxl>0 for nvar values of l.
11:   xldx* – Real (Kind=nag_wp) arrayInput
Note: the second dimension of the array x must be at least max1,m.
On entry: if mode='S' the kth row of x must contain xk. That is xkl must contain the kth sample value for the lth variable, for k=1,2,,nobs and l=1,2,,m. Otherwise x is not referenced.
12:   ldx – IntegerInput
On entry: the first dimension of the array x as declared in the (sub)program from which g03dbf is called.
Constraints:
  • if mode='S', ldxnobs;
  • otherwise ldx1.
13:   dlddng – Real (Kind=nag_wp) arrayOutput
On exit: the squared distances.
If mode='S', dkj contains the squared distance of the kth sample point from the jth group mean, Dkj2, for k=1,2,,nobs and j=1,2,,ng.
If mode='M' and equal='U', dij contains the squared distance between the ith mean and the jth mean, Dij2, for i=1,2,,ng and j=1,2,,i-1,i+1,,ng. The elements dii are not referenced, for i=1,2,,ng.
If mode='M' and equal='E', dij contains the squared distance between the ith mean and the jth mean, Dij2, for i=1,2,,ng and j=1,2,,i-1. Since Dij=Dji the elements dij are not referenced, for i=1,2,,ng and j=i+1,,ng.
14:   ldd – IntegerInput
On entry: the first dimension of the array d as declared in the (sub)program from which g03dbf is called.
Constraints:
  • if mode='S', lddnobs;
  • if mode='M', lddng.
15:   wk2×nvar – Real (Kind=nag_wp) arrayWorkspace
16:   ifail – IntegerInput/Output
On entry: ifail must be set to 0, -1 or 1. If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation for details.
For environments where it might be inappropriate to halt program execution when an error is detected, the value -1 or 1 is recommended. If the output of error messages is undesirable, then the value 1 is recommended. Otherwise, if you are not familiar with this argument, the recommended value is 0. When the value -1 or 1 is used it is essential to test the value of ifail on exit.
On exit: ifail=0 unless the routine detects an error or a warning has been flagged (see Section 6).

6
Error Indicators and Warnings

If on entry ifail=0 or -1, explanatory error messages are output on the current error message unit (as defined by x04aaf).
Errors or warnings detected by the routine:
ifail=1
On entry, equal=value.
Constraint: equal='E' or 'U'.
On entry, ldd=value and ng=value.
Constraint: lddng.
On entry, ldd=value and nobs=value.
Constraint: lddnobs.
On entry, ldgmn=value and ng=value.
Constraint: ldgmnng.
On entry, ldx=value and nobs=value.
Constraint: ldxnobs.
On entry, m=value and nvar=value.
Constraint: mnvar.
On entry, mode=value.
Constraint: mode='M' or 'S'.
On entry, ng=value.
Constraint: ng2.
On entry, nobs=value.
Constraint: nobs1.
On entry, nvar=value.
Constraint: nvar1.
ifail=2
On entry, diagonal element value of R=0.
On entry, diagonal element value of Rj=0 for j=value.
On entry, nvar=value and value values of isx>0.
Constraint: exactly nvar elements of isx>0.
ifail=-99
An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.
ifail=-399
Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.
ifail=-999
Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

7
Accuracy

The accuracy will depend upon the accuracy of the input R or Rj matrices.

8
Parallelism and Performance

g03dbf makes calls to BLAS and/or LAPACK routines, which may be threaded within the vendor library used by this implementation. Consult the documentation for the vendor library for further information.
Please consult the X06 Chapter Introduction for information on how to control and interrogate the OpenMP environment used within this routine. Please also consult the Users' Note for your implementation for any additional implementation-specific information.

9
Further Comments

If the distances are to be used for discrimination, see also g03dcf.

10
Example

The data, taken from Aitchison and Dunsmore (1975), is concerned with the diagnosis of three ‘types’ of Cushing's syndrome. The variables are the logarithms of the urinary excretion rates (mg/24hr) of two steroid metabolites. Observations for a total of 21 patients are input and the group means and R matrices are computed by g03daf. A further six observations of unknown type are input, and the distances from the group means of the 21 patients of known type are computed under the assumption that the within-group variance-covariance matrices are not equal. These results are printed and indicate that the first four are close to one of the groups while observations 5 and 6 are some distance from any group.

10.1
Program Text

Program Text (g03dbfe.f90)

10.2
Program Data

Program Data (g03dbfe.d)

10.3
Program Results

Program Results (g03dbfe.r)