Module_Mul_Stat_Procedures

Copyright 2024 IRD

This file is part of statpack.

statpack is free software: you can redistribute it and/or modify it under the terms of the GNU Lesser General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.

statpack is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more details.

You can find a copy of the GNU Lesser General Public License in the statpack/doc directory.

Authors: Pascal Terray (LOCEAN/IPSL, Paris, France)


MODULE EXPORTING SUBROUTINES AND FUNCTIONS FOR MULTIVARIATE STATISTICAL COMPUTATIONS

LATEST REVISION : 12/03/2024


subroutine comp_cor ( x, y, first, last, xstat, ystat, xycor, xyn, z, prob, ndf_max, cortest, cov )

Purpose

COMP_COR computes estimates of mean, variance and correlation coefficient from two data vectors XX and YY.

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, input subvector containing size(X) observations from the vector of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data vector XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X) observations from another vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = size( X ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current subvector is the first subvector of the data vector XX (or YY).
  • FIRST = false the current subvector is not the first subvector of the data vector XX (or YY).
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current subvector is the last subvector of the data vector XX (or YY).
  • LAST = false the current subvector is not the last subvector of the data vector XX (or YY).
XSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. XSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, XSTAT contains the following statistics:

  • XSTAT(1) contains the mean value of the data vector XX.
  • XSTAT(2) contains the variance of the data vector XX.

The size of XSTAT must verify: size( XSTAT ) = 2.

YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. YSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR. XYCOR should not be changed between calls to COMP_COR.

On exit, when LAST=true, XYCOR contains the correlation coefficient between XX(:) and YY(:).

XYN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYN contains count of observations from previous calls to COMP_COR. XYN should not be changed between calls to COMP_COR.

On exit, XYN contains the number of observations in the data vectors XX and YY.

Z (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

PROB (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, PROB gives the probability that the random sample of XYN observation pairs YY(:) and XX(:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB.
  • If XYN-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB.

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

CORTEST (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, a probability. CORTEST is the sum of the areas (equal) in both tails of the Student’s t distribution with XYN-2 degrees of freedom.

CORTEST must verify: 0. < P < 1.

On exit, the two-tail CORTEST quantile of the sample correlation coefficient, that is a value R such that the probability of the absolute value of a sample correlation coefficient computed from XYN observation pairs being greater than R is CORTEST under the null hypothesis of no correlation in the parent normal population.

CORTEST needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficient between the data vectors XX and YY is computed instead of correlation coefficient. If COV=true, Z and PROB are set to Nan code.

COV needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficient with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to Nan code.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor ( x, y, first, last, xstat, ystat, xycor, xyn, dimvar, z, prob, ndf_max, cortest, cov )

Purpose

COMP_COR computes estimates of means, variances and correlation coefficients (with a data vector YY) from a data matrix XX.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which basic univariate and bivariate statistics are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3-DIMVAR) observations from a vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3-DIMVAR ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix XX.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix XX.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix XX.
  • LAST = false the current submatrix is not the last submatrix of the data matrix XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. XSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,1) contains the mean values of the data matrix XX.
  • XSTAT(:,2) contains the variances of the data matrix XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ),
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. YSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR. XYCOR should not be changed between calls to COMP_COR.

On exit, when LAST=true, XYCOR(i) contains the correlation coefficient between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:).

The size of XYCOR must verify: size( XYCOR ) = size( X, DIMVAR ).

XYN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYN contains count of observations from previous calls to COMP_COR. XYN should not be changed between calls to COMP_COR.

On exit, XYN contains the number of observations in the data matrix XX and the data vector YY.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The size of Z must verify: size( Z ) = size( X, DIMVAR ).

PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, PROB(i) gives the probability that the random sample of XYN observation pairs YY(:) and XX(i,:) ( XX(:,i) if DIMVAR=2 ) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i).
  • If XYN-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

CORTEST (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, a probability. CORTEST is the sum of the areas (equal) in both tails of the Student’s t distribution with XYN-2 degrees of freedom.

CORTEST must verify: 0. < P < 1.

On exit, the two-tail CORTEST quantile of the sample correlation coefficient, that is a value R such that the probability of the absolute value of a sample correlation coefficient computed from XYN observation pairs being greater than R is CORTEST under the null hypothesis of no correlation in the parent normal population.

CORTEST needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)
On entry, if COV=true, covariance coefficients between the data matrix XX and data vector YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to Nan code. COV needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to Nan code.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor ( x, y, first, last, xstat, ystat, xycor, xyn, z, prob, ndf_max, cortest, cov )

Purpose

COMP_COR computes estimates of means, variances and correlation coefficients (with a data vector YY) from a data tridimensional array XX.

Arguments

X (INPUT) real(stnd), dimension(:,:,:)
On entry, input subarray containing size(X,3) observations on size(X,1) by size(X,2) variables from the tridimensional array of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data array XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3) observations from a vector of data YY for which correlation coefficients with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3 ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current subarray is the first subarray of the data array XX.
  • FIRST = false the current subarray is not the first subarray of the data array XX.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current subarray is the last subarray of the data array XX.
  • LAST = false the current subarray is not the last subarray of the data array XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. XSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,:,1) contains the mean values of the data array XX.
  • XSTAT(:,:,2) contains the variances of the data array XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, 1 ) ,
  • size( XSTAT, 2 ) = size( X, 2 ) ,
  • size( XSTAT, 3 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. YSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR. XYCOR should not be changed between calls to COMP_COR.

On exit, when LAST=true, XYCOR(i,j) contains the correlation coefficient between XX(i,j,:) and YY(:).

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, 1 ) ,
  • size( XYCOR, 2 ) = size( X, 2 ) .
XYN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYN contains count of observations from previous calls to COMP_COR. XYN should not be changed between calls to COMP_COR.

On exit, XYN contains the number of observations in the data array XX and the data vector YY.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The shape of Z must verify:

  • size( Z, 1 ) = size( X, 1 ) ,
  • size( Z, 2 ) = size( X, 2 ) .
PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, PROB(i,j) gives the probability that the random sample of XYN observation pairs YY(:) and XX(i,j,:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The shape of PROB must verify:

  • size( PROB, 1 ) = size( X, 1 ) ,
  • size( PROB, 2 ) = size( X, 2 ) .
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i,j).
  • If XYN-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i,j).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

CORTEST (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, a probability. CORTEST is the sum of the areas (equal) in both tails of the Student’s t distribution with XYN-2 degrees of freedom.

CORTEST must verify: 0. < P < 1.

On exit, the two-tail CORTEST quantile of the sample correlation coefficient, that is a value R such that the probability of the absolute value of a sample correlation coefficient computed from XYN observation pairs being greater than R is CORTEST under the null hypothesis of no correlation in the parent normal population.

CORTEST needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)
On entry, if COV=true, covariance coefficients between the data matrices XX and YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to Nan code. COV needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to Nan code.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor ( x, y, first, last, xstat, ystat, xycor, xyn, dimvar, dimvary, z, prob, ndf_max, cortest, cov )

Purpose

COMP_COR computes estimates of means, variances and correlation coefficients between two data matrices YY and XX.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which basic univariate and bivariate statistics are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the matrix of data YY for which basic univariate and bivariate statistics are desired. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVAR ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. XSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the data matrix XX.
  • XSTAT(:,2) contains the variances of the data matrix XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR. YSTAT should not be changed between calls to COMP_COR.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the data matrix YY.
  • YSTAT(:,2) contains the variances of the data matrix YY.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 2 .
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR. XYCOR should not be changed between calls to COMP_COR.

On exit, when LAST=true, XYCOR(i,j) contains the correlation coefficient between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVAR=2 and DIMVARY=2 ).

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) .
XYN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR (e.g., when FIRST=true), XYN contains count of observations from previous calls to COMP_COR. XYN should not be changed between calls to COMP_COR.

On exit, XYN contains the number of observations in the data matrices XX and YY.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVAR = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The shape of Z must verify:

  • size( Z, 1 ) = size( X, DIMVAR ) ,
  • size( Z, 2 ) = size( Y, DIMVARY ) .
PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, PROB(i,j) gives the probability that the random sample of XYN observation pairs XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVAR=2 and DIMVARY=2 ) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

The shape of PROB must verify:

  • size( PROB, 1 ) = size( X, DIMVAR ) ,
  • size( PROB, 2 ) = size( Y, DIMVARY ) .
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i,j).
  • If XYN-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i,j).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

CORTEST (INPUT/OUTPUT, OPTIONAL) real(stnd)

On entry, a probability. CORTEST is the sum of the areas (equal) in both tails of the Student’s t distribution with XYN-2 degrees of freedom.

CORTEST must verify: 0. < P < 1.

On exit, the two-tail CORTEST quantile of the sample correlation coefficient, that is a value R such that the probability of the absolute value of a sample correlation coefficient computed from XYN observation pairs being greater than R is CORTEST under the null hypothesis of no correlation in the parent normal population.

CORTEST needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficients between the data matrices XX and YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to Nan code.

COV needs to be specified only on the last call to COMP_COR (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to Nan code.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss ( x, y, first, last, xstat, ystat, xycor, xymiss, z, prob, ndf_max, cov )

Purpose

COMP_COR_MISS computes estimates of mean, variance and correlation coefficient from two data vectors XX and YY possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, input subvector containing size(X) observations from the vector of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data vector XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X) observations from another vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = size( X ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current subvector is the first subvector of the data vector XX (or YY).
  • FIRST = false the current subvector is not the first subvector of the data vector XX (or YY).
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current subvector is the last subvector of the data vector XX (or YY).
  • LAST = false the current subvector is not the last subvector of the data vector XX (or YY).
XSTAT (INPUT/OUTPUT) real(stnd), dimension(4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XSTAT contains the following statistics:

  • XSTAT(1) contains the mean value of the data vector XX.
  • XSTAT(2) contains the variance of the data vector XX.
  • XSTAT(3) contains the the number of non-missing observations in the data vector XX.
  • XSTAT(4) is used as workspace.

The size of XSTAT must verify: size( XSTAT ) = 4.

YSTAT (INPUT/OUTPUT) real(stnd), dimension(4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. YSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.
  • YSTAT(3) contains the the number of non-missing observations in the data vector YY.
  • YSTAT(4) is used as workspace.

The size of YSTAT must verify: size( YSTAT ) = 4.

XYCOR (INPUT/OUTPUT) real(stnd), dimension(4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XYCOR should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(1) contains the correlation coefficient between XX(:) and YY(:).
  • XYCOR(2) contains the incidence value between XX(:) and YY(:). XYCOR(2) indicates the number of non-missing pairs of observations which were used in the calculation of XYCOR(1).
  • XYCOR(3:4) is used as workspace.

The size of XYCOR must verify: size( XYCOR ) = 4.

XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics are computed on all the observations where XX and YY are not missing (see Further Details).
Z (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR(1).

Z needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

PROB (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, PROB gives the probability that the random sample of XYCOR(2) observation pairs YY(:) and XX(:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows, if:

  • XYCOR(2)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB.
  • XYCOR(2)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB.

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficient between the data vectors XX and YY are computed instead of correlation coefficient. If COV=true, Z and PROB are set to XYMISS .

COV needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficient with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to XYMISS.

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss ( x, y, first, last, xstat, ystat, xycor, xymiss, dimvar, z, prob, ndf_max, cov )

Purpose

COMP_COR_MISS computes estimates of means, variances and correlation coefficients (with another data vector YY) from a data matrix XX possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which basic univariate and bivariate statistics are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3-DIMVAR) observations from a vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3-DIMVAR ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix XX.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix XX.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix XX.
  • LAST = false the current submatrix is not the last submatrix of the data matrix XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,1) contains the mean values of the data matrix XX.
  • XSTAT(:,2) contains the variances of the data matrix XX.
  • XSTAT(:,3) contains the the numbers of non-missing observations in the data matrix XX.
  • XSTAT(:,4) is used as workspace.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 4 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. YSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.
  • YSTAT(3) contains the the number of non-missing observations in the data vector YY.
  • YSTAT(4) is used as workspace.

The size of YSTAT must verify: size( YSTAT ) = 4.

XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XYCOR should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(i,1) contains the correlation coefficients between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:).
  • XYCOR(i,2) contains the incidence values between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:). XYCOR(i,2) indicates the numbers of non-missing pairs of observations which were used in the calculation of XYCOR(i,1).
  • XYCOR(:,3:4) is used as workspace.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XYCOR, 2 ) = 4 .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics are computed on all the observations where X and Y are not missing (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR(:,1).

Z needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The size of Z must verify: size( Z ) = size( X, DIMVAR ).

PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, PROB(i) gives the probability that the random sample of XYCOR(i,2) observation pairs YY(:) and XX(i,:) ( XX(:,i) if DIMVAR=2 ) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYCOR(i,2)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i).
  • If XYCOR(i,2)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficients between the data matrix XX and data vector YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to XYMISS .

COV needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the statistics are set to XYMISS .

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss ( x, y, first, last, xstat, ystat, xycor, xymiss, z, prob, ndf_max, cov )

Purpose

COMP_COR_MISS computes estimates of means, variances and correlation coefficients (with another data vector YY) from a data tridimensional array XX possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:,:)
On entry, input subarray containing size(X,3) observations on size(X,1) by size(X,2) variables from the array of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data array XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3) observations from a vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3 ) .

FIRST (INPUT) logical(lgl)

On entry, if

  • FIRST = true the current subarray is the first subarray of the data array XX.
  • FIRST = false the current subarray is not the first subarray of the data array XX.
LAST (INPUT) logical(lgl)

On entry, if

  • LAST = true the current subarray is the last subarray of the data array XX.
  • LAST = false the current subarray is not the last subarray of the data array XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,:,1) contains the mean values of the data array XX.
  • XSTAT(:,:,2) contains the variances of the data array XX.
  • XSTAT(:,:,3) contains the numbers of non-missing observations in the data array XX.
  • XSTAT(:,:,4) is used as workspace.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, 1 ) ,
  • size( XSTAT, 2 ) = size( X, 2 ) ,
  • size( XSTAT, 3 ) = 4 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. YSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.
  • YSTAT(3) contains the the number of non-missing observations in the data vector YY.
  • YSTAT(4) is used as workspace.

The size of YSTAT must verify: size( YSTAT ) = 4.

XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XYCOR should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(i,j,1) contains the correlation coefficients between XX(i,j,:) and YY(:).
  • XYCOR(i,j,2) contains the incidence values between XX(i,j,:) and YY(:). XYCOR(i,j,2) indicates the numbers of valid pairs of observations which were used in the calculation of XYCOR(:,:,1).
  • XYCOR(:,:,3:4) is used as workspace.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, 1 ) ,
  • size( XYCOR, 2 ) = size( X, 2 ) ,
  • size( XYCOR, 3 ) = 4 .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics are computed on all the observations where X and Y are not missing (see Further Details).
Z (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR(:,:,1).

Z needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The shape of Z must verify:

  • size( Z, 1 ) = size( X, 1 ) ,
  • size( Z, 2 ) = size( X, 2 ) .
PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, PROB(i,j) gives the probability that the random sample of XYCOR(i,j,2) observation pairs YY(:) and XX(i,j,:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The shape of PROB must verify:

  • size( PROB, 1 ) = size( X, 1 ) ,
  • size( PROB, 2 ) = size( X, 2 ) .
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYCOR(i,j,2)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i,j).
  • If XYCOR(i,j,2)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i,j).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficients between the data array XX and the data vector YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to XYMISS .

COV needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the statistics are set to XYMISS .

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss ( x, y, first, last, xstat, ystat, xycor, xymiss, dimvar, dimvary, z, prob, ndf_max, cov )

Purpose

COMP_COR_MISS computes estimates of means, variances and correlation coefficients between two data matrices YY and XX possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which basic univariate and bivariate statistics are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the matrix of data YY for which basic univariate and bivariate statistics are desired. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVAR ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the data matrix XX.
  • XSTAT(:,2) contains the variances of the data matrix XX.
  • XSTAT(:,3) contains the the numbers of non-missing observations in the data matrix XX.
  • XSTAT(:,4) is used as workspace.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 4 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. YSTAT should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the data matrix YY.
  • YSTAT(:,2) contains the variances of the data matrix YY.
  • YSTAT(:,3) contains the the numbers of non-missing observations in the data matrix YY.
  • YSTAT(:,4) is used as workspace.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 4 .
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:,4)

On entry, after the first call to COMP_COR_MISS (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS. XYCOR should not be changed between calls to COMP_COR_MISS.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(i,j,1) contains the correlation coefficients between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVAR=2 and DIMVARY=2 ).
  • XYCOR(i,j,2) contains the incidence values between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVAR=2 and DIMVARY=2). XYCOR(i,j,2) indicates the numbers of non-missing pairs of observations which were used in the calculation of XYCOR(i,j,1).
  • XYCOR(:,:,3:4) is used as workspace.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) ,
  • size( XYCOR, 3 ) = 4 .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics are computed on all the observations where X and Y are not missing (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVAR = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The shape of Z must verify:

  • size( Z, 1 ) = size( X, DIMVAR ) ,
  • size( Z, 2 ) = size( Y, DIMVARY ) .
PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, PROB(i,j) gives the probability that the random sample of XYCOR(i,j,2) observation pairs XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVAR=2 and DIMVARY=2 ) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

The shape of PROB must verify:

  • size( PROB, 1 ) = size( X, DIMVAR ) ,
  • size( PROB, 2 ) = size( Y, DIMVARY ) .
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYCOR(i,j,2)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i,j).
  • If XYCOR(i,j,2)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i,j).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV=true, covariance coefficients between the data matrices XX and YY are computed instead of correlation coefficients. If COV=true, Z and PROB are set to XYMISS .

COV needs to be specified only on the last call to COMP_COR_MISS (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the statistics are set to XYMISS .

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation and regression analysis, see

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss2 ( x, y, first, last, xstat, ystat, xycor, xyn, xymiss, z, prob, ndf_max )

Purpose

COMP_COR_MISS2 computes estimates of mean, variance and correlation coefficient from two data vectors XX and YY possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:)
On entry, input subvector containing size(X) observations from the vector of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data vector XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X) observations from another vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = size( X ) .

FIRST (INPUT) logical(lgl)

On entry, if

  • FIRST = true the current subvector is the first subvector of the data vector XX (or YY).
  • FIRST = false the current subvector is not the first subvector of the data vector XX (or YY).
LAST (INPUT) logical(lgl)

On entry, if

  • LAST = true the current subvector is the last subvector of the data vector XX (or YY).
  • LAST = false the current subvector is not the last subvector of the data vector XX (or YY).
XSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XSTAT should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XSTAT contains the following statistics:

  • XSTAT(1) contains the mean value of the data vector XX.
  • XSTAT(2) contains the variance of the data vector XX.

The size of XSTAT must verify: size( XSTAT ) = 2.

YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. YSTAT should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, YSTAT contains the following statistics:

  • YSTAT(1) contains the mean value of the data vector YY.
  • YSTAT(2) contains the variance of the data vector YY.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XYCOR should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XYCOR contains the correlation coefficient between XX(:) and YY(:).

XYN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYN contains count of valid pairs of observations from previous calls to COMP_COR_MISS2. XYN should not be changed between calls to COMP_COR_MISS2.

On exit, XYN contains the incidence value between XX(:) and YY(:). XYN indicates the number of non-missing pairs of observations which were used in the calculation of XSTAT, YSTAT and XYCOR.

XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics are computed on all valid pairs of observations where XX and YY are not missing (see Further Details).
Z (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR.

Z needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

PROB (OUTPUT, OPTIONAL) real(stnd)

On exit, when LAST=true, PROB gives the probability that the random sample of XYN observation pairs YY(:) and XX(:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB.
  • If XYN-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB.

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficient with only one pass through the data.

If fewer than two valid observations were present, the pertinent statistics are set to XYMISS.

The univariate and bivariate statistics are computed from all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss2 ( x, y, first, last, xstat, ystat, xycor, xyn, xymiss, dimvar, z, prob, ndf_max )

Purpose

COMP_COR_MISS2 computes estimates of means, variances and correlation coefficients (with another data vector YY) from a data matrix XX possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which basic univariate and bivariate statistics are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3-DIMVAR) observations from a vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3-DIMVAR ) .

FIRST (INPUT) logical(lgl)

On entry, if

  • FIRST = true the current submatrix is the first submatrix of the data matrix XX.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix XX.
LAST (INPUT) logical(lgl)

On entry, if

  • LAST = true the current submatrix is the last submatrix of the data matrix XX.
  • LAST = false the current submatrix is not the last submatrix of the data matrix XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XSTAT should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,1) contains the mean values of the data matrix XX.
  • XSTAT(:,2) contains the variances of the data matrix XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 2.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. YSTAT should not be changed between calls to COMP_COR_MISS2.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( YSTAT, 2 ) = 2.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XYCOR should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XYCOR(i) contains the correlation coefficient between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:).

The size of XYCOR must verify: size( XYCOR ) = size( X, DIMVAR ).

XYN (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYN contains counts of valid pairs of observations from previous calls to COMP_COR_MISS2. XYN should not be changed between calls to COMP_COR_MISS2.

On exit, XYN(i) contains the incidence value between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:). XYN(i) indicates the number of non-missing pairs of observations which were used in the calculation of XSTAT and XYCOR.

The size of XYN must verify: size( XYN ) = size( X, DIMVAR ).

XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics for variable XX(i,:) ( XX(:,i) if DIMVAR=2 ) are computed on all valid pairs of observations where XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:) are not missing (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

Z (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR(:).

Z needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

The size of Z must verify: size( Z ) = size( X, DIMVAR ).

PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true, PROB(i) gives the probability that the random sample of XYN(i) observation pairs YY(:) and XX(i,:) ( XX(:,i) if DIMVAR=2 ) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN(i)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i).
  • If XYN(i)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the statistics are set to XYMISS .

The univariate and bivariate statistics are computed from all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cor_miss2 ( x, y, first, last, xstat, ystat, xycor, xyn, xymiss, z, prob, ndf_max )

Purpose

COMP_COR_MISS2 computes estimates of means, variances and correlation coefficients (with another data vector YY) from a data tridimensional array XX possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:,:)
On entry, input subarray containing size(X,3) observations on size(X,1) by size(X,2) variables from the array of data XX for which basic univariate and bivariate statistics are desired. If all the data are available at once, X can be the full data array XX.
Y (INPUT) real(stnd), dimension(:)

On entry, input subvector containing size(X,3) observations from a vector of data YY for which correlation coefficient with XX is desired. If all the data are available at once, YY can be the full data vector.

The size of Y must verify: size( Y ) = SIZE( X, 3 ) .

FIRST (INPUT) logical(lgl)

On entry, if

  • FIRST = true the current subarray is the first subarray of the data array XX.
  • FIRST = false the current subarray is not the first subarray of the data array XX.
LAST (INPUT) logical(lgl)

On entry, if

  • LAST = true the current subarray is the last subarray of the data array XX.
  • LAST = false the current subarray is not the last subarray of the data array XX.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XSTAT should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XSTAT contains the following statistics on all variables:

  • XSTAT(:,:,1) contains the mean values of the data array XX.
  • XSTAT(:,:,2) contains the variances of the data array XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, 1 ) ,
  • size( XSTAT, 2 ) = size( X, 2 ) ,
  • size( XSTAT, 3 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. YSTAT should not be changed between calls to COMP_COR_MISS2.

The shape of XSTAT must verify:

  • size( YSTAT, 1 ) = size( X, 1 ) ,
  • size( YSTAT, 2 ) = size( X, 2 ) ,
  • size( YSTAT, 3 ) = 2.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_COR_MISS2. XYCOR should not be changed between calls to COMP_COR_MISS2.

On exit, when LAST=true, XYCOR(i,j) contains the correlation coefficient between XX(i,j,:) and YY(:).

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, 1 ) ,
  • size( XYCOR, 2 ) = size( X, 2 ) .
XYN (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_COR_MISS2 (e.g., when FIRST=true), XYN contains counts of valid pairs of observations from previous calls to COMP_COR_MISS2. XYN should not be changed between calls to COMP_COR_MISS2. On exit, XYN(i,j) contains the incidence value between XX(i,j,:) and YY(:). XYN(i,j) indicates the number of non-missing pairs of observations which were used in the calculation of XSTAT and XYCOR.

The shape of XYN must verify:

  • size( XYN, 1 ) = size( X, 1 ) ,
  • size( XYN, 2 ) = size( X, 2 ) .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate and bivariate statistics for variable XX(i,j,:) are computed on all valid pairs of observations where XX(i,j,:) and YY(:) are not missing (see Further Details).
Z (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, Z contains the Fisher’s Z transformation of XYCOR(:,:).

Z needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

The shape of Z must verify:

  • size( Z, 1 ) = size( X, 1 ) ,
  • size( Z, 2 ) = size( X, 2 ) .
PROB (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, PROB(i,j) gives the probability that the random sample of XYN(i,j) observation pairs YY(:) and XX(i,j,:) came from a bivariate normal population with a correlation coefficient equal to zero.

PROB needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

The shape of PROB must verify:

  • size( PROB, 1 ) = size( X, 1 ) ,
  • size( PROB, 2 ) = size( X, 2 ) .
NDF_MAX (INPUT, OPTIONAL) integer(i4b)

On entry, when argument PROB is present, NDF_MAX is used as follows:

  • If XYN(i,j)-2 is lower or equal to NDF_MAX, the t_density is integrated for computing PROB(i,j).
  • If XYN(i,j)-2 is greater than NDF_MAX, an asymptotic series is used for computing PROB(i,j).

The default is 20.

NDF_MAX needs to be specified only on the last call to COMP_COR_MISS2 (e.g., when LAST=true).

Further Details

The subroutine computes the basic univariate statistics and the correlation coefficients with only one pass through the data.

If fewer than two valid observations were present, the statistics are set to XYMISS .

The univariate and bivariate statistics are computed from all valid pairs of observations.

For more details on correlation and regression analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine update_cor ( xstat, ystat, xycor, xyn, xstat2, ystat2, xycor2, xyn2 )

Purpose

UPDATE_COR computes sample means and corrected sums of squares and cross-products for a sample of size XYN+XYN2 given the means and corrected sum of squares and cross-products for two subsamples of size XYN and XYN2 as output by a call to COMP_COR when LAST=false on the two subsamples separetely.

The sample means, variances and coefficient correlation for the sample of size XYN+XYN2 may be obtained by a call to COMP_COR with LAST=true.

Arguments

XSTAT (INPUT/OUTPUT) real(stnd), dimension(2)
On entry, the XSTAT argument of COMP_COR for the first subsample. On exit, the XSTAT argument of the combined sample.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the first subsample. On exit, the YSTAT argument of the combined sample.
XYCOR (INPUT/OUTPUT) real(stnd)
On entry, the XYCOR argument of COMP_COR for the first subsample. On exit, the XYCOR argument of the combined sample.
XYN (INPUT/OUTPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the first subsample. On exit, the XYN argument of the combined sample.
XSTAT2 (INPUT) real(stnd), dimension(2)
On entry, the XSTAT argument of COMP_COR for the second subsample.
YSTAT2 (INPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the second subsample.
XYCOR2 (INPUT) real(stnd)
On entry, the XYCOR argument of COMP_COR for the second subsample.
XYN2 (INPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the second subsample.

Further Details

One possible application of this subroutine is to parallel processing. If one has two or more processors available, the sample can be split up into smaller subsamples, and the means and corrected sums of squares and cross-products computed for each subsample independently using COMP_COR. The means and corrected sums of squares and cross-products for the original sample can then be calculated using UPDATE_COR. The means, variances and correlation coefficient for the original sample can be computed by a final call to COMP_COR with LAST=true.

This subroutine is adapted from:

  1. Chan, T.F., Golub, G.H, and Leveque, R.J., 1979:
    Updating formulae and a pairwise algorithm for computing sample variances. STAN-CS-79-773.

subroutine update_cor ( xstat, ystat, xycor, xyn, xstat2, ystat2, xycor2, xyn2 )

Purpose

UPDATE_COR computes sample means and corrected sums of squares and cross-products for a sample of size XYN+XYN2 given the means and corrected sums of squares and cross-products for two subsamples of size XYN and XYN2 as output by a call to COMP_COR when LAST=false on the two subsamples separetely.

The sample means, variances and coefficient correlations for the sample of size XYN+XYN2 may be obtained by a call to COMP_COR with LAST=true.

Arguments

XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)
On entry, the XSTAT argument of COMP_COR for the first subsample. On exit, the XSTAT argument of the combined sample.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the first subsample. On exit, the YSTAT argument of the combined sample.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the XYCOR argument of COMP_COR for the first subsample. On exit, the XYCOR argument of the combined sample.
XYN (INPUT/OUTPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the first subsample. On exit, the XYN argument of the combined sample.
XSTAT2 (INPUT) real(stnd), dimension(:,2)
On entry, the XSTAT argument of COMP_COR for the second subsample.
YSTAT2 (INPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the second subsample.
XYCOR2 (INPUT) real(stnd), dimension(:)
On entry, the XYCOR argument of COMP_COR for the second subsample.
XYN2 (INPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the second subsample.

Further Details

One possible application of this subroutine is to parallel processing. If one has two or more processors available, the sample can be split up into smaller subsamples, and the means and corrected sums of squares and cross-products computed for each subsample independently using COMP_COR. The means and corrected sums of squares and cross-products for the original sample can then be calculated using UPDATE_COR. The means, variances and correlation coefficient for the original sample can be computed by a final call to COMP_COR with LAST=true.

This subroutine is adapted from:

  1. Chan, T.F., Golub, G.H, and Leveque, R.J., 1979:
    Updating formulae and a pairwise algorithm for computing sample variances. STAN-CS-79-773.

subroutine update_cor ( xstat, ystat, xycor, xyn, xstat2, ystat2, xycor2, xyn2 )

Purpose

UPDATE_COR computes sample means and corrected sums of squares and cross-products for a sample of size XYN+XYN2 given the means and corrected sums of squares and cross-products for two subsamples of size XYN and XYN2 as output by a call to COMP_COR when LAST=false on the two subsamples separetely.

The sample means, variances and coefficient correlations for the sample of size XYN+XYN2 may be obtained by a call to COMP_COR with LAST=true.

Arguments

XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)
On entry, the XSTAT argument of COMP_COR for the first subsample. On exit, the XSTAT argument of the combined sample.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the first subsample. On exit, the YSTAT argument of the combined sample.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the XYCOR argument of COMP_COR for the first subsample. On exit, the XYCOR argument of the combined sample.
XYN (INPUT/OUTPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the first subsample. On exit, the XYN argument of the combined sample.
XSTAT2 (INPUT) real(stnd), dimension(:,:,2)
On entry, the XSTAT argument of COMP_COR for the second subsample.
YSTAT2 (INPUT) real(stnd), dimension(2)
On entry, the YSTAT argument of COMP_COR for the second subsample.
XYCOR2 (INPUT) real(stnd), dimension(:,:)
On entry, the XYCOR argument of COMP_COR for the second subsample.
XYN2 (INPUT) real(stnd)
On entry, the XYN argument of COMP_COR for the second subsample.

Further Details

One possible application of this subroutine is to parallel processing. If one has two or more processors available, the sample can be split up into smaller subsamples, and the means and corrected sums of squares and cross-products computed for each subsample independently using COMP_COR. The means and corrected sums of squares and cross-products for the original sample can then be calculated using UPDATE_COR. The means, variances and correlation coefficient for the original sample can be computed by a final call to COMP_COR with LAST=true.

This subroutine is adapted from:

  1. Chan, T.F., Golub, G.H, and Leveque, R.J., 1979:
    Updating formulae and a pairwise algorithm for computing sample variances. STAN-CS-79-773.

subroutine update_cor_miss2 ( xstat, ystat, xycor, xyn, xstat2, ystat2, xycor2, xyn2 )

Purpose

UPDATE_COR_MISS2 computes sample means and corrected sums of squares and cross-products for a sample of size XYN+XYN2, possibly containing missing values, given the means and corrected sums of squares and cross-products for two subsamples of size XYN and XYN2 as output by a call to COMP_COR_MISS2 when LAST=false on the two subsamples separetely.

The sample means, variances and coefficient correlations for the sample of size XYN+XYN2 may be obtained by a call to COMP_COR_MISS2 with LAST=true.

Arguments

XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)
On entry, the XSTAT argument of COMP_COR_MISS2 for the first subsample. On exit, the XSTAT argument of the combined sample.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)
On entry, the YSTAT argument of COMP_COR_MISS2 for the first subsample. On exit, the YSTAT argument of the combined sample.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the XYCOR argument of COMP_COR_MISS2 for the first subsample. On exit, the XYCOR argument of the combined sample.
XYN (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the XYN argument of COMP_COR_MISS2 for the first subsample. On exit, the XYN argument of the combined sample.
XSTAT2 (INPUT) real(stnd), dimension(:,2)
On entry, the XSTAT argument of COMP_COR_MISS2 for the second subsample.
YSTAT2 (INPUT) real(stnd), dimension(:,2)
On entry, the YSTAT argument of COMP_COR_MISS2 for the second subsample.
XYCOR2 (INPUT) real(stnd), dimension(:)
On entry, the XYCOR argument of COMP_COR_MISS2 for the second subsample.
XYN2 (INPUT) real(stnd), dimension(:)
On entry, the XYN argument of COMP_COR_MISS2 for the second subsample.

Further Details

One possible application of this subroutine is to parallel processing. If one has two or more processors available, the sample can be split up into smaller subsamples, and the means and corrected sums of squares and cross-products computed for each subsample independently using COMP_COR_MISS2. The means and corrected sums of squares and cross-products for the original sample can then be calculated using UPDATE_COR_MISS2. The means, variances and correlation coefficient for the original sample can be computed by a final call to COMP_COR_MISS2 with LAST=true.

This subroutine is adapted from:

  1. Chan, T.F., Golub, G.H, and Leveque, R.J., 1979:
    Updating formulae and a pairwise algorithm for computing sample variances. STAN-CS-79-773.

subroutine update_cor_miss2 ( xstat, ystat, xycor, xyn, xstat2, ystat2, xycor2, xyn2 )

Purpose

UPDATE_COR_MISS2 computes sample means and corrected sums of squares and cross-products for a sample of size XYN+XYN2, possibly containing missing values, given the means and corrected sums of squares and cross-products for two subsamples of size XYN and XYN2 as output by a call to COMP_COR_MISS2 when LAST=false on the two subsamples separetely.

The sample means, variances and coefficient correlations for the sample of size XYN+XYN2 may be obtained by a call to COMP_COR_MISS2 with LAST=true.

Arguments

XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)
On entry, the XSTAT argument of COMP_COR_MISS2 for the first subsample. On exit, the XSTAT argument of the combined sample.
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,:,2)
On entry, the YSTAT argument of COMP_COR_MISS2 for the first subsample. On exit, the YSTAT argument of the combined sample.
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the XYCOR argument of COMP_COR_MISS2 for the first subsample. On exit, the XYCOR argument of the combined sample.
XYN (INPUT/OUTPUT) real(stnd), dimension(:,:)
On entry, the XYN argument of COMP_COR_MISS2 for the first subsample. On exit, the XYN argument of the combined sample.
XSTAT2 (INPUT) real(stnd), dimension(:,:,2)
On entry, the XSTAT argument of COMP_COR_MISS2 for the second subsample.
YSTAT2 (INPUT) real(stnd), dimension(:,:,2)
On entry, the YSTAT argument of COMP_COR_MISS2 for the second subsample.
XYCOR2 (INPUT) real(stnd), dimension(:,:)
On entry, the XYCOR argument of COMP_COR_MISS2 for the second subsample.
XYN2 (INPUT) real(stnd), dimension(:,:)
On entry, the XYN argument of COMP_COR_MISS2 for the second subsample.

Further Details

One possible application of this subroutine is to parallel processing. If one has two or more processors available, the sample can be split up into smaller subsamples, and the means and corrected sums of squares and cross-products computed for each subsample independently using COMP_COR_MISS2. The means and corrected sums of squares and cross-products for the original sample can then be calculated using UPDATE_COR_MISS2. The means, variances and correlation coefficient for the original sample can be computed by a final call to COMP_COR_MISS2 with LAST=true.

This subroutine is adapted from:

  1. Chan, T.F., Golub, G.H, and Leveque, R.J., 1979:
    Updating formulae and a pairwise algorithm for computing sample variances. STAN-CS-79-773.

subroutine permute_cor ( x, y, xstat, ystat, xycor, prob, nrep, initseed )

Purpose

PERMUTE_COR performs a permutation test of a correlation coefficient between two data vectors Y and X.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the input data vector X. On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.
Y (INPUT) real(stnd), dimension(:)

On entry, the input data vector Y.

The size of Y must verify: SIZE( Y ) = SIZE( X ) .

XSTAT (INPUT) real(stnd), dimension(2)

On entry, XSTAT must contain the following statistics as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(1) contains the mean value of the data vector X.
  • XSTAT(2) contains the variance of the data vector X.

The size of XSTAT must verify: size( XSTAT ) = 2.

YSTAT (INPUT) real(stnd), dimension(2)

On entry, YSTAT must contain the following statistics, as output by COMP_COR subroutine in argument YSTAT:

  • YSTAT(1) contains the mean value of the data vector Y.
  • YSTAT(2) contains the variance of the data vector Y.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT) real(stnd)
On entry, XYCOR contains the correlation coefficient between X(:) and Y(:). XYCOR must be specified as output by COMP_COR subroutine.
PROB (OUTPUT) real(stnd)
On exit, PROB gives the critical probability associated with XYCOR, as computed by a permutation test with NREP shuffles.
NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the permutation test of the correlation coefficient.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details

This subroutine is parallelized if OPENMP is used.

For more details and algorithm, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapters 5 and 8, 484 pp., ISBN:9780521012300
  2. Noreen, E.W., 1989:
    Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3

subroutine permute_cor ( x, y, xstat, ystat, xycor, prob, dimvar, nrep, initseed )

Purpose

PERMUTE_COR performs permutation tests of correlation coefficients between a data vector Y and a data matrix X.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, input matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables for which permutation tests are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details.

On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.

Y (INPUT) real(stnd), dimension(:)

On entry, input vector containing size(X,3-DIMVAR) observations for which permutation tests are desired.

The size of Y must verify: size( Y ) = SIZE( X, 3-DIMVAR ) .

XSTAT (INPUT) real(stnd), dimension(:,2)

On entry, XSTAT must contain the following statistics on all variables, as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(:,1) contains the mean values of the data matrix X.
  • XSTAT(:,2) contains the variances of the data matrix X.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT) real(stnd), dimension(2)

On entry, YSTAT must contain the following statistics, as output by COMP_COR subroutine in argument YSTAT:

  • YSTAT(1) contains the mean value of the data vector Y.
  • YSTAT(2) contains the variance of the data vector Y.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT) real(stnd), dimension(:)

On entry, XYCOR(i) contains the correlation coefficient between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:). XYCOR must be specified as output by COMP_COR subroutine.

The size of XYCOR must verify: size( XYCOR ) = size( X, DIMVAR ).

PROB (OUTPUT) real(stnd), dimension(:)

On exit, PROB(i) gives the critical probability associated with XYCOR(i), as computed by a permutation test with NREP shuffles.

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows:

  • DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input matrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the permutation test of the correlation coefficients.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details

This subroutine is parallelized if OPENMP is used.

For more details and algorithm, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapters 5 and 8, 484 pp., ISBN:9780521012300
  2. Noreen, E.W., 1989:
    Computer-intensive methods for testing hypotheses: an introduction. Wiley and Sons, New York, USA, ISBN:978-0-471-61136-3

subroutine phase_scramble_cor ( x, y, xstat, ystat, xycor, prob, nrep, method, norm, initseed )

Purpose

PHASE_SCRAMBLE_COR performs phase-scrambled bootstrap tests of a correlation coefficient between two data vectors Y and X.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:)
On entry, the input data vector X. On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.
Y (INPUT) real(stnd), dimension(:)

On entry, the input data vector Y.

The size of Y must verify: SIZE( Y ) = SIZE( X ) .

XSTAT (INPUT) real(stnd), dimension(2)

On entry, XSTAT must contain the following statistics as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(1) contains the mean value of the data vector X.
  • XSTAT(2) contains the variance of the data vector X.

The size of XSTAT must verify: size( XSTAT ) = 2.

YSTAT (INPUT) real(stnd), dimension(2)

On entry, YSTAT must contain the following statistics, as output by COMP_COR subroutine in argument YSTAT:

  • YSTAT(1) contains the mean value of the data vector Y.
  • YSTAT(2) contains the variance of the data vector Y.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT) real(stnd)
On entry, XYCOR contains the correlation coefficient between X(:) and Y(:). XYCOR must be specified as output by COMP_COR subroutine.
PROB (OUTPUT) real(stnd)
On exit, PROB gives the critical probability associated with XYCOR, as computed by a phase-scrambled bootstrap test with NREP shuffles.
NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the phase-scrambled bootstrap test of the correlation coefficient.

The default is 99.

METHOD (INPUT, OPTIONAL) integer(i4b)

On entry, determine the phase randomisation algorithm used to generate surrogate series.

On entry, if

  • METHOD = 1 : the phase randomisation algorithm of Theiler is used;
  • METHOD = 2 : the phase randomisation algorithm of Davison and Hinkley is used.

The default is METHOD = 1.

NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if NORM=true, then normal margins are used in the phase-scrambled algorithm, otherwise exact empirical margins are used.

The default is NORM=true.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details

This subroutine is parallelized if OPENMP is used.

The tests are adapted from:

  1. Ebisuzaki, W., 1997:
    A method to estimate the statistical significance of a correlation when the data are serially correlated. Journal of climate, vol. 10, 2147-2153.
  2. Davison, A.C., and Hinkley, D.V., 1997:
    Bootstrap methods and their application. Cambridge University Press, Cambridge, UK. doi:10.1017/CBO9780511802843
  3. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, J.D., 1992:
    Testing for nonlinearity in time series: the method of surrogate data. Physica D, vol. 58, 77-94, doi:10.1016/0167-2789(92)90102-s
  4. Braun, W.J., and Kulperger, R.J., 1997:
    Properties of a fourier bootstrap method for time series. Communications in Statistics - Theory and Methods, vol 26, 1329-1336, doi:10.1080/03610929708831985

subroutine phase_scramble_cor ( x, y, xstat, ystat, xycor, prob, dimvar, nrep, method, norm, initseed )

Purpose

PHASE_SCRAMBLE_COR performs phase-scrambled bootstrap tests of correlation coefficients between a data vector Y and a data matrix X.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, input matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables for which phase-scrambled bootstrap tests are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details.

On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.

Y (INPUT) real(stnd), dimension(:)

On entry, input vector containing size(X,3-DIMVAR) observations for which phase-scrambled bootstrap tests are desired.

The size of Y must verify: size( Y ) = SIZE( X, 3-DIMVAR ) .

XSTAT (INPUT) real(stnd), dimension(:,2)

On entry, XSTAT must contain the following statistics on all variables, as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(:,1) contains the mean values of the data matrix X.
  • XSTAT(:,2) contains the variances of the data matrix X.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT) real(stnd), dimension(2)

On entry, YSTAT must contain the following statistics, as output by COMP_COR subroutine in argument YSTAT:

  • YSTAT(1) contains the mean value of the data vector Y.
  • YSTAT(2) contains the variance of the data vector Y.

The size of YSTAT must verify: size( YSTAT ) = 2.

XYCOR (INPUT) real(stnd), dimension(:)

On entry, XYCOR(i) contains the correlation coefficient between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:). XYCOR must be specified as output by COMP_COR subroutine.

The size of XYCOR must verify: size( XYCOR ) = size( X, DIMVAR ).

PROB (OUTPUT) real(stnd), dimension(:)

On exit, PROB(i) gives the critical probability associated with XYCOR(i), as computed by a phase-scrambled bootstrap test with NREP shuffles.

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input matrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the phase-scrambled bootstrap test of the correlation coefficients.

The default is 99.

METHOD (INPUT, OPTIONAL) integer(i4b)

On entry, determine the phase randomisation algorithm used to generate surrogate series.

On entry, if

  • METHOD = 1 : the phase randomisation algorithm of Theiler is used;
  • METHOD = 2 : the phase randomisation algorithm of Davison and Hinkley is used.

The default is METHOD = 1.

NORM (INPUT, OPTIONAL) logical(lgl)

On entry, if NORM=true, then normal margins are used in the phase-scrambled algorithm, otherwise exact empirical margins are used.

The default is NORM=true.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

Further Details

This subroutine is parallelized if OPENMP is used.

The tests are adapted from:

  1. Ebisuzaki, W., 1997:
    A method to estimate the statistical significance of a correlation when the data are serially correlated. Journal of climate, vol. 10, 2147-2153.
  2. Davison, A.C., and Hinkley, D.V., 1997:
    Bootstrap methods and their application. Cambridge University Press, Cambridge, UK. doi:10.1017/CBO9780511802843
  3. Theiler, J., Eubank, S., Longtin, A., Galdrikian, B., and Farmer, J.D., 1992:
    Testing for nonlinearity in time series: the method of surrogate data. Physica D, vol. 58, 77-94, doi:10.1016/0167-2789(92)90102-s
  4. Braun, W.J., and Kulperger, R.J., 1997:
    Properties of a fourier bootstrap method for time series. Communications in Statistics - Theory and Methods, vol 26, 1329-1336, doi:10.1080/03610929708831985

subroutine bootstrap_cor ( x, y, xstat, xycor, prob, nrep, initseed, periodicity, season, block_size )

Purpose

BOOTSTRAP_COR performs a moving block bootstrap test of a correlation coefficient between two data vectors X and Y.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, the input data vector X.

On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.

Y (INPUT) real(stnd), dimension(:)

On entry, the input data vector Y.

The size of Y must verify: SIZE( Y ) = SIZE( X ) .

XSTAT (INPUT) real(stnd), dimension(2)

On entry, XSTAT must contain the following statistics as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(1) contains the mean value of the data vector X.
  • XSTAT(2) contains the variance of the data vector X.

The size of XSTAT must verify: size( XSTAT ) = 2.

XYCOR (INPUT) real(stnd)
On entry, XYCOR contains the correlation coefficient between X(:) and Y(:). XYCOR must be specified as output by COMP_COR subroutine.
PROB (OUTPUT) real(stnd)
On exit, PROB gives the critical probability associated with XYCOR, as computed by a moving block bootstrap test with NREP shuffles.
NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the moving block bootstrap test of the correlation coefficient.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

PERIODICITY (INPUT, OPTIONAL) integer(i4b)

On entry, argument PERIODICITY specifies that the indice, i, of the first observation of each selected block in the moving block bootstrap algorithm verifies the condition i=1+(PERIODICITY * j) where j is a random positive integer. PERIODICITY must be greater than zero and less than size(X).

By default, PERIODICITY is set to 1.

SEASON (INPUT, OPTIONAL) integer(i4b)

On entry, argument SEASON specifies that the input time series is a repetition of the same season for different years and SEASON specifies the length of the season. SEASON must be greater than zero and size(X) must be a multiple of SEASON. If the optional argument PERIODICITY is used, SEASON must also be greater or equal to PERIODICITY.

By default, SEASON is set to size(X).

BLOCK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, argument BLOCK_SIZE specifies the size of the block in the moving block bootstrap. BLOCK_SIZE must be greater than zero and less than size(X). If the optional argument PERIODICITY is used, BLOCK_SIZE must also be greater or equal to PERIODICITY. Moreover, if the optional argument SEASON is used, BLOCK_SIZE must also be less than SEASON.

By default, BLOCK_SIZE is set to 1 or to PERIODICITY if this optional argument is used.

Further Details

This subroutine is parallelized if OPENMP is used.

The test is adapted from:

  1. Davison, A.C., and Hinkley, D.V., 1997:
    Bootstrap methods and their application. Cambridge University Press, Cambridge, UK. doi:10.1017/CBO9780511802843

subroutine bootstrap_cor ( x, y, xstat, xycor, prob, dimvar, nrep, initseed, periodicity, season, block_size )

Purpose

BOOTSTRAP_COR performs a moving block bootstrap test of a correlation coefficients between a data vector Y and a data matrix X.

Arguments

X (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, input matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables for which moving block bootstrap tests are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details.

On exit, the data are standardized with the univariate statistics stored in the XSTAT argument.

Y (INPUT) real(stnd), dimension(:)

On entry, input vector containing size(X,3-DIMVAR) observations for which moving block bootstrap tests are desired.

The size of Y must verify: size( Y ) = size( X, 3-DIMVAR ) .

XSTAT (INPUT) real(stnd), dimension(:,2)

On entry, XSTAT must contain the following statistics on all variables, as output by COMP_COR subroutine in argument XSTAT:

  • XSTAT(:,1) contains the mean values of the data matrix X.
  • XSTAT(:,2) contains the variances of the data matrix X.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVAR ) ,
  • size( XSTAT, 2 ) = 2 .
XYCOR (INPUT) real(stnd), dimension(:)

On entry, XYCOR(i) contains the correlation coefficient between XX(i,:) ( XX(:,i) if DIMVAR=2 ) and YY(:). XYCOR must be specified as output by COMP_COR subroutine.

The size of XYCOR must verify: size( XYCOR ) = size( X, DIMVAR ).

PROB (OUTPUT) real(stnd), dimension(:)

On exit, PROB(i) gives the critical probability associated with XYCOR(i), as computed by a moving block bootstrap test with NREP shuffles.

The size of PROB must verify: size( PROB ) = size( X, DIMVAR ).

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input matrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

NREP (INPUT, OPTIONAL) integer(i4b)

On entry, when argument NREP is present, NREP specifies the number of shuffles for the moving block bootstrap test of the correlation coefficient.

The default is 99.

INITSEED (INPUT, OPTIONAL) logical(lgl)

On entry, if INITSEED=true, a call to RANDOM_SEED_() without arguments is done in the subroutine, in order to initiates a non-repeatable reset of the seed used by the STATPACK random generator.

The default is INITSEED=false.

PERIODICITY (INPUT, OPTIONAL) integer(i4b)

On entry, argument PERIODICITY specifies that the indice, i, of the first observation of each selected block in the moving block bootstrap algorithm verifies the condition i=1+(PERIODICITY * j) where j is a random positive integer. PERIODICITY must be greater than zero and less than size(X,3-DIMVAR).

By default, PERIODICITY is set to 1.

SEASON (INPUT, OPTIONAL) integer(i4b)

On entry, argument SEASON specifies that the input time series is a repetition of the same season for different years and SEASON specifies the length of the season. SEASON must be greater than zero and size(X) must be a multiple of SEASON. If the optional argument PERIODICITY is used, SEASON must also be greater or equal to PERIODICITY.

By default, SEASON is set to size(X,3-DIMVAR).

BLOCK_SIZE (INPUT, OPTIONAL) integer(i4b)

On entry, argument BLOCK_SIZE specifies the size of the block in the moving block bootstrap. BLOCK_SIZE must be greater than zero and less than size(X). If the optional argument PERIODICITY is used, BLOCK_SIZE must also be greater or equal to PERIODICITY. Moreover, if the optional argument SEASON is used, BLOCK_SIZE must also be less than SEASON.

By default, BLOCK_SIZE is set to 1 or to PERIODICITY if this optional argument is used.

Further Details

This subroutine is parallelized if OPENMP is used.

The test is adapted from:

  1. Davison, A.C., and Hinkley, D.V., 1997:
    Bootstrap methods and their application. Cambridge University Press, Cambridge, UK. doi:10.1017/CBO9780511802843

subroutine comp_cormat ( x, first, last, xmean, xcor, xn, dimvar, xstd, cov, fill, failure )

Purpose

COMP_CORMAT computes estimates of means and variance-covariance or correlation matrix from a data matrix.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which means, variances and covariances are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XMEAN (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), XMEAN contains the variable means from previous calls to COMP_CORMAT. XMEAN should not be changed between calls to COMP_CORMAT.

On exit, when LAST=true, XMEAN contains the variable means

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), the matrix XCOR contains the upper triangle of the corrected sums of squared and cross-products matrix computed from previous calls to COMP_CORMAT. XCOR should not be changed between calls to COMP_CORMAT.

On exit, when LAST=true, XCOR contains the upper triangle of the symetric correlation or variance-covariance matrix as controlled by the COV argument. If the optional argument FILL is present and equal to true, the lower triangle of XCOR is also filled.

The shape of XCOR must verify:

  • size( XCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XCOR, 2 ) = size( X, DIMVAR ) .
XN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), XN contains count of observations from previous calls to COMP_CORMAT. XN should not be changed between calls to COMP_CORMAT.

On exit, XN contains the number of observations in the data matrix XX.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSTD needs to be specified only on the last call to COMP_CORMAT (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XCOR contains the variances-covariances matrix, when LAST=true.
  • COV= false, XCOR contains the correlation matrix, when LAST=true.

By default, the correlation matrix is output.

COV needs to be specified only on the last call to COMP_CORMAT (e.g., when LAST=true).

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows, if:

  • FILL= true, the lower triangle of XCOR is filled, when LAST=true.
  • FILL= false, the lower triangle of XCOR is not filled, when LAST=true.

By default, the lower triangle of XCOR is not filled.

FILL needs to be specified only on the last call to COMP_CORMAT (e.g., when LAST=true).

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit, when argument FAILURE is present:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that the observations on some variable were constant and the correlations were requested.

Further Details

The subroutine computes the means and correlation matrix with only one pass through the data.

If the observations on some variable were constant, the pertinent correlations are set to nan() code .

If fewer than two valid observations were present, the correlations are set to nan() code .

If fewer than one valid observation is present, the means are also set to nan() code .

For more details on correlation analysis, see

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cormat ( x, first, last, xmean, xcorp, xn, dimvar, xstd, cov, failure )

Purpose

COMP_CORMAT computes estimates of means and variance-covariance or correlation matrix from a data matrix.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which means, variances and covariances are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XMEAN (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), XMEAN contains the variable means from previous calls to COMP_CORMAT. XMEAN should not be changed between calls to COMP_CORMAT. On exit, when LAST=true, XMEAN contains the variable means

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XCORP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), the linear array XCORP contains the upper triangle of the corrected sums of squared and cross-products matrix, packed columnwise, computed from previous calls to COMP_CORMAT. XCORP should not be changed between calls to COMP_CORMAT.

On exit, when LAST=true, XCORP contains the correlation or variance-covariance matrix as controlled by the COV argument. XCORP is stored in symmetric storage mode (see further details).

The size of XCORP must verify: size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2.

XN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_CORMAT (e.g., when when FIRST=true), XN contains count of observations from previous calls to COMP_CORMAT. XN should not be changed between calls to COMP_CORMAT.

On exit, XN contains the number of observations in the data matrix XX.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSTD needs to be specified only on the last call to COMP_CORMAT (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XCORP contains the variances-covariances matrix, when LAST=true.
  • COV= false, XCORP contains the correlation matrix, when LAST=true.

By default, the correlation matrix is output.

COV needs to be specified only on the last call to COMP_CORMAT (e.g., when LAST=true).

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit, when argument FAILURE is present:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that the observations on some variable were constant and the correlations were requested.

Further Details

The subroutine computes the means and the correlation matrix with only one pass through the data.

On exit, the upper triangle of the symmetric correlation or variance-covariance matrix XCOR is packed columnwise in the linear array XCORP. More precisely, the j-th column of XCOR is stored in the array XCORP as follows:

XCORP(i + (j-1) * j/2) = XCOR(i,j) for 1<=i<=j;

If the observations on some variable were constant, the pertinent correlations are set to nan() code .

If fewer than two valid observations were present, the correlations are set to nan() code .

If fewer than one valid observation is present, the means are also set to nan() code .

For more details on correlation analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cormat_miss ( x, first, last, xmean, xcor, xn, xmiss, dimvar, xstd, cov, fill, failure )

Purpose

COMP_CORMAT_MISS computes estimates of means and variance-covariance or correlation matrix from a data matrix possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which means, variances and covariances are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XMEAN (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), XMEAN(:,1) contains the variable means from previous calls to COMP_CORMAT_MISS. XMEAN should not be changed between calls to COMP_CORMAT_MISS.

On exit, when LAST=true, XMEAN(:,1) contains the variable means computed from all non-missing observations in the data matrix. XMEAN(:,2) is used as workspace.

The shape of XMEAN must verify:

  • size( XMEAN, 1 ) = size( X, DIMVAR ) ,
  • size( MEAN, 2 ) = 2 .
XCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), the matrix XCOR contains the upper triangle of the corrected sums of squared and cross-products matrix computed from previous calls to COMP_CORMAT_MISS. XCOR should not be changed between calls to COMP_CORMAT_MISS.

On exit, when LAST=true, XCOR contains the upper triangle of the symetric correlation or variance-covariance matrix as controlled by the COV argument. If the optional argument FILL is present and equal to true, the lower triangle of XCOR is also filled.

The shape of XCOR must verify:

  • size( XCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XCOR, 2 ) = size( X, DIMVAR ) .
XN (INPUT/OUTPUT) real(stnd), dimension(:,3)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), XN is used as workspace to accumulate quantities from previous calls to COMP_CORMAT_MISS. XN should not be changed between calls to COMP_CORMAT_MISS.

On exit, XN(:,1) contains the upper triangle of the matrix of the incidence values between each pair of variables, packed columnwise, in a linear array. XN(i + (j-1) * j/2,1) indicates the numbers of non-missing pairs which were used in the calculation of XCOR(i,j) for 1<=i<=j . XN(:,2:3) is used as workspace.

The shape of XN must verify:

  • size( XN, 1 ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2 ,
  • size( XN, 2 ) = 3 .
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The means and the correlations are computed on all the observations where X are not missing (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSTD needs to be specified only on the last call to COMP_CORMAT_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XCORP contains the variances-covariances matrix, when LAST=true.
  • COV= false, XCORP contains the correlation matrix, when LAST=true.

By default, the correlation matrix is output.

COV needs to be specified only on the last call to COMP_CORMAT_MISS (e.g., when LAST=true).

FILL (INPUT, OPTIONAL) logical(lgl)

On entry, when argument FILL is present, FILL is used as follows, if:

  • FILL= true, the lower triangle of XCOR is filled, when LAST=true.
  • FILL= false, the lower triangle of XCOR is not filled, when LAST=true.

By default, the lower triangle of XCOR is not filled.

FILL needs to be specified only on the last call to COMP_CORMAT_MISS (e.g., when LAST=true).

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit, when argument FAILURE is present:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some pair of variables or that the observations on some variable were constant and the correlations were requested.

Further Details

The subroutine computes the means and the correlation matrix with only one pass through the data.

If the observations on some variable were constant, the pertinent correlations are set to XMISS.

If fewer than two valid observations were present for some pair of variables, the pertinent correlations are set to XMISS.

If fewer than one valid observation is present for some variables, the pertinent means are also set to XMISS.

The means and standard-deviations of the data matrix are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_cormat_miss ( x, first, last, xmean, xcorp, xn, xmiss, dimvar, xstd, cov, failure )

Purpose

COMP_CORMAT_MISS computes estimates of means and variance-covariance or correlation matrix from a data matrix possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which means, variances and covariances are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XMEAN (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), XMEAN(:,1) contains the variable means from previous calls to COMP_CORMAT_MISS. XMEAN should not be changed between calls to COMP_CORMAT_MISS.

On exit, when LAST=true, XMEAN(:,1) contains the variable means computed from all non-missing observations in the data matrix. XMEAN(:,2) is used as workspace.

The shape of XMEAN must verify:

  • size( XMEAN, 1 ) = size( X, DIMVAR ) ,
  • size( MEAN, 2 ) = 2 .
XCORP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), the linear array XCORP contains the upper triangle of the corrected sums of squared and cross-products matrix, packed columnwise, computed from previous calls to COMP_CORMAT_MISS. XCORP should not be changed between calls to COMP_CORMAT_MISS.

On exit, when LAST=true, XCORP contains the correlation or variance-covariance matrix as controlled by the COV argument. XCORP is stored in symmetric storage mode (see further details).

The size of XCORP must verify: size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2.

XN (INPUT/OUTPUT) real(stnd), dimension(:,3)

On entry, after the first call to COMP_CORMAT_MISS (e.g., when FIRST=true), XN is used as workspace to accumulate quantities from previous calls to COMP_CORMAT_MISS. XN should not be changed between calls to COMP_CORMAT_MISS.

On exit, XN(:,1) contains the incidence values between each pair of variables. XN(i,1) indicates the numbers of non-missing pairs of observations which were used in the calculation of XCORP(i).

The shape of XN must verify:

  • size( XN, 1 ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2 ,
  • size( XN, 2 ) = 3 .
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The means and the correlations are computed on all the observations where X are not missing (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSTD needs to be specified only on the last call to COMP_CORMAT_MISS (e.g., when LAST=true).

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XCORP contains the variances-covariances matrix, when LAST=true.
  • COV= false, XCORP contains the correlation matrix, when LAST=true.

By default, the correlation matrix is output.

COV needs to be specified only on the last call to COMP_CORMAT_MISS (e.g., when LAST=true).

FAILURE (OUTPUT, OPTIONAL) logical(lgl)

On exit, when argument FAILURE is present:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some variable or that the observations on some variable were constant and the correlations were requested.

Further Details

The subroutine computes the means and the correlation matrix with only one pass through the data.

On exit, the upper triangle of the symmetric correlation or variance-covariance matrix XCOR is packed columnwise in the linear array XCORP. More precisely, the j-th column of XCOR is stored in the array XCORP as follows:

XCORP(i + (j-1) * j/2) = XCOR(i,j) for 1<=i<=j;

If the observations on some variable were constant, the pertinent correlations are set to XMISS .

If fewer than two valid observations were present on some variable, the pertinent correlations are set to XMISS.

If fewer than one valid observation is present for some variable, the pertinent means are also set to XMISS.

The means and standard-deviations of the data matrix are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations.

For more details on correlation analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 8, 484 pp., ISBN:9780521012300

subroutine comp_eof ( x, first, last, xeigval, xeigvec, xn, failure, dimvar, cov, sort, maxiter, xmean, xstd, xeigvar, xcorp )

Purpose

COMP_EOF computes estimates of Empirical Orthogonal Functions (EOF; also known as Principal Component Analysis) from a data matrix.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Empirical Orthogonal Functions are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XEIGVAL (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_EOF (e.g., when FIRST=true), XEIGVAL contains temporary results from previous calls to COMP_EOF. XEIGVAL should not be changed between calls to COMP_EOF.

On exit, when LAST=true, XEIGVAL contains the eigenvalues of th variance-covariance (or correlation) matrix from the data matrix. The near zero eigenvalues made negative by round off errors are set to zero.

The size of XEIGVAL must verify: size( XEIGVAL ) = size( X, DIMVAR ).

XEIGVEC (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_EOF (e.g., when FIRST=true), the matrix XEIGVEC contains temporary results from previous calls to COMP_EOF. XEIGVEC should not be changed between calls to COMP_EOF.

On exit, when LAST=true, XEIGVEC contains the eigenvectors of th variance-covariance (or correlation) matrix from the data matrix.

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) = size( X, DIMVAR ) .
XN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_EOF (e.g., when FIRST=true), XN contains count of observations from previous calls to COMP_EOF. XN should not be changed between calls to COMP_EOF.

On exit, XN contains the number of observations in the data matrix.

FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the eigenvectors and the eigenvalues.

On exit when LAST=false, FAILURE is always set to false.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XEIGVEC contains the eigenvectors of the variances-covariances matrix, when LAST=true.
  • COV= false, XEIGVEC contains the eigenvectors of the correlation matrix, when LAST=true.

By default, the eigenvectors of the correlation matrix are output.

COV needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

SORT needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of the covariance matrix . The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(XEIGVAL). Convergence usually occurs in about 2 * size(XEIGVAL) QR sweeps.

The default is 30.

MAXITER needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

XMEAN (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XMEAN is present, XMEAN contains the variable means.

XMEAN needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

XSTD needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XEIGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XEIGVAR is present, XEIGVAR contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL.

XEIGVAR needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

The size of XEIGVAR must verify: size( XEIGVAR ) = size( X, DIMVAR ).

XCORP (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XCORP is present, XCORP contains the upper triangle of the correlation or variance-covariance matrix, as controlled by the COV argument, stored in symmetric storage mode. The upper triangle of the symmetric correlation or variance-covariance matrix is packed columnwise in the linear array XCORP. More precisely, the j-th column of this matrix is stored in the array XCORP as follows:

XCORP(i + (j-1) * j/2,2) = XCOR(i,j) for 1<=i<=j;

XCORP needs to be specified only on the last call to COMP_EOF (e.g., when LAST=true).

The size of XCORP must verify:

  • size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2

Further Details

The subroutine computes the Empirical Orthogonal Functions with only one pass through the data.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVAR) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XEIGVAL, XEIGVEC, XEIGVAR and XCORP are globally set to NaN code.

For more details on EOF or PCA analysis, see:

  1. Jackson, J.E., 2003:
    A user’s guide to principal components. John Wiley and Sons, New York, USA, 592 pp., ISBN:978-0-471-47134-9
  2. Jolliffe, I.T., 2002:
    Principal component analysis. Springer-Verlag, New York, USA, 487 pp., 2nd Ed, ISBN:978-0-387-22440-4
  3. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_eof2 ( x, first, last, xeigval, xcorp, xn, failure, dimvar, cov, savecor, maxiter, ortho, xmean, xstd, xeigvar, xeigvec )

Purpose

COMP_EOF2 computes estimates of Empirical Orthogonal Functions (EOF; also known as Principal Component Analysis) from a data matrix.

COMP_EOF2 computes all the eigenvalues, and optionally selected eigenvectors (by inverse iteration), of the covariance (or correlation matrix) from a data matrix.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Empirical Orthogonal Functions are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XEIGVAL (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_EOF2 (e.g., when FIRST=true), XEIGVAL contains temporary results from previous calls to COMP_EOF2. XEIGVAL should not be changed between calls to COMP_EOF2.

On exit, when LAST=true, XEIGVAL contains the eigenvalues of th variance-covariance (or correlation) matrix from the data matrix. The near zero eigenvalues made negative by round off errors are set to zero. The eigenvalues are sorted in descending order. The eigenvectors in XEIGVEC are reordered accordingly.

The size of XEIGVAL must verify: size( XEIGVAL ) = size( X, DIMVAR ).

XCORP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_EOF2 (e.g., when FIRST=true), the linear array XCORP contains the upper triangle of the corrected sums of squared and cross-products matrix, packed columnwise, computed from previous calls to COMP_EOF2. XCORP should not be changed between calls to COMP_EOF2.

On exit, when LAST=true and SAVECOR=true, XCORP contains the correlation or variance-covariance matrix as controlled by the COV argument. XCORP is stored in symmetric storage mode (see further details). If SAVECOR=false, the correlation matrix is not saved on exit. In this case XCORP does not contain useful information.

The size of XCORP must verify: size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2.

XN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_EOF2 (e.g., when FIRST=true), XN contains count of observations from previous calls to COMP_EOF2. XN should not be changed between calls to COMP_EOF2.

On exit, XN contains the number of observations in the data matrix.

FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the eigenvalues or that some eigenvectors failed to converge with MAXITER inverse iterations.

On exit when LAST=false, FAILURE is always set to false.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, the eigenvalues and the eigenvectors are computed from the variances-covariances matrix, when LAST=true.
  • COV= false, the eigenvalues and the eigenvectors are computed from the correlation matrix, when LAST=true.

By default, the eigenvalues and eigenvectors of the correlation matrix are computed.

COV needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

SAVECOR (INPUT, OPTIONAL) logical(lgl)

On exit, when argument SAVECOR is present and LAST=true, SAVECOR is used as follows, if:

  • SAVECOR= true, the correlation (or covariance) matrix is saved in packed form in argument XCORP.
  • SAVECOR= false, the correlation (or covariance) matrix is destroyed.

By default, the correlation (or covariance) matrix is destroyed.

SAVECOR needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine for computing the eigenvectors. By default, 2 inverse iterations are performed for all the eigenvectors. This optional argument is used only if the optional argument XEIGVEC is present.

MAXITER needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed eigenvectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the optional argument XEIGVEC is present.

The default is FALSE.

ORTHO needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

XMEAN (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XMEAN is present, XMEAN contains the variable means.

XMEAN needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ) .

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

XSTD needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XEIGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XEIGVAR is present, XEIGVAR contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL.

XEIGVAR needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

The size of XEIGVAR must verify: size( XEIGVAR ) = size( X, DIMVAR ).

XEIGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, XEIGVEC contains the first size(XEIGVEC,2) eigenvectors of the variance-covariance (or correlation) matrix from the data matrix.

XEIGVEC needs to be specified only on the last call to COMP_EOF2 (e.g., when LAST=true).

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) <= size( X, DIMVAR ) .

Further Details

The subroutine computes the means and the covariance (or correlation) matrix with only one pass through the data.

On exit, if SAVECOR= true, the upper triangle of the symmetric correlation or variance-covariance matrix XCOR is packed columnwise in the linear array XCORP. More precisely, the j-th column of XCOR is stored in the array XCORP as follows:

XCORP(i + (j-1) * j/2) = XCOR(i,j) for 1<=i<=j;

Eigenvalues and selected eigenvectors are computed from the packed correlation matrix when LAST=true.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVAR) = 0), in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XEIGVAL, XEIGVEC, XEIGVAR and XCORP are globally set to NaN code.

For more details on EOF or PCA analysis, see:

  1. Jackson, J.E., 2003:
    A user’s guide to principal components. John Wiley and Sons, New York, USA, 592 pp., ISBN:978-0-471-47134-9
  2. Jolliffe, I.T., 2002:
    Principal component analysis. Springer-Verlag, New York, USA, 487 pp., 2nd Ed, ISBN:978-0-387-22440-4
  3. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_eof3 ( x, dimvar, failure, xcorp, xeigval, xeigvec, maxiter, ortho )

Purpose

COMP_EOF3 computes estimates of Empirical Orthogonal Functions (EOF; also known as Principal Component Analysis) from a data matrix X with n observations.

COMP_EOF3 computes the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X, all the eigenvalues, and selected eigenvectors (by inverse iteration), of this matrix product.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input data matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables for which Empirical Orthogonal Functions or Principal Components are desired.
DIMVAR (INPUT) integer(i4b)

On entry, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables and the matrix product (1/n) (X * X’) is computed.
  • DIMVAR = 2, the input matrix X contains size(X,1) observations on size(X,2) variables and the matrix product (1/n) (X’ * X) is computed.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that maximum accuracy was not achieved when computing the eigenvalues or that some eigenvectors failed to converge with MAXITER inverse iterations.
XCORP (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, XCORP(:) contains the matrix product (1/n) (X’ * X) or (1/n) (X * X’), stored in symmetric storage mode.

The upper triangle of the symmetric matrix product matrix is packed columnwise in the linear array XCORP(:). More precisely, the j-th column of this matrix is stored in the array XCORP(:) as follows:

XCORP(i + (j-1) * j/2,1) = XCOR(i,j) for 1<=i<=j;

The size of XCORP must verify: size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2

XEIGVAL (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit:

  • XEIGVAL(:,1) contains the eigenvalues in decreasing order of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X. The near zero eigenvalues made negative by round off errors are set to zero.
  • XEIGVAL(:,2) contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL(:,1).

The shape of XEIGVAL must verify:

  • size( XEIGVAL, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVAL, 2 ) = 2 .
XEIGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, XEIGVEC contains the first size(XEIGVEC,2) eigenvectors of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X.

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) <= size( X, DIMVAR ) .
MAXITER (INPUT, OPTIONAL) integer(i4b)
The number of inverse iterations performed in the subroutine for computing the eigenvectors. By default, 2 inverse iterations are performed for all the eigenvectors. This optional argument is used only if the XEIGVEC is present.
ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed eigenvectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the XEIGVEC is present.

The default is FALSE.

Further Details

The subroutine computes the Empirical Orthogonal Functions or the Principal Components with only one pass through the data.

If size(X,3-DIMVAR)<=0, the subroutine set FAILURE to true and returns without doing any computations.

For more details on EOF or PCA analysis, see:

  1. Jackson, J.E., 2003:
    A user’s guide to principal components. John Wiley and Sons, New York, USA, 592 pp., ISBN:978-0-471-47134-9
  2. Jolliffe, I.T., 2002:
    Principal component analysis. Springer-Verlag, New York, USA, 487 pp., 2nd Ed, ISBN:978-0-387-22440-4
  3. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_eof_miss ( x, first, last, xeigval, xeigvec, xcorp, xmiss, failure, dimvar, cov, sort, maxiter, xmean, xstd )

Purpose

COMP_EOF_MISS computes estimates of Empirical Orthogonal Functions (EOF; also known as Principal Component Analysis) from a data matrix possibly containing missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Empirical Orthogonal Functions are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XEIGVAL (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_EOF_MISS (e.g., when FIRST=true), XEIGVAL contains temporary results from previous calls to COMP_EOF_MISS. XEIGVAL should not be changed between calls to COMP_EOF_MISS. On exit, when LAST=true:

  • XEIGVAL(:,1) contains the eigenvalues of the variance-covariance (or correlation) matrix from the data matrix. The near zero eigenvalues made negative by round off errors or because the variance-covariance (or correlation) matrix from the data matrix with missing values is not positive definite are set to zero.
  • XEIGVAL(:,2) contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL(:,1).

The shape of XEIGVAL must verify:

  • size( XEIGVAL, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVAL, 2 ) = 2 .
XEIGVEC (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_EOF_MISS (e.g., when FIRST=true), the matrix XEIGVEC contains temporary results from previous calls to COMP_EOF_MISS. XEIGVEC should not be changed between calls to COMP_EOF_MISS.

On exit, when LAST=true, XEIGVEC contains the eigenvectors of th variance-covariance (or correlation) matrix from the data matrix.

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) = size( X, DIMVAR ) .
XCORP (INPUT/OUTPUT) real(stnd), dimension(:,3)

On entry, after the first call to COMP_EOF_MISS (e.g., when FIRST=true), XCORP is used as workspace to accumulate quantities from previous calls to COMP_EOF_MISS. XCORP should not be changed between calls to COMP_EOF_MISS.

On exit, when LAST=true:

  • XCORP(:,1) contains the correlation or variance-covariance matrix, as controlled by the COV argument, stored in symmetric storage mode. The upper triangle of the symmetric correlation or variance-covariance matrix is packed columnwise in the linear array XCORP(:,1). More precisely, the j-th column of this matrix is stored in the array XCORP(:,1) as follows:

    XCORP(i + (j-1) * j/2,1) = XCOR(i,j) for 1<=i<=j;

  • XCORP(:,2) contains the upper triangle of the matrix of the incidence values between each pair of variables, packed columnwise, in a linear array. XCORP(i + (j-1) * j/2,2) indicates the numbers of non-missing pairs which were used in the calculation of the covariance (or correlation) between variables i and j, for 1<=i<=j .

  • XCORP(:,3) is used as workspace and contains no useful information.

The shape of XCORP must verify:

  • size( XCORP, 1 ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2 ,
  • size( XCORP, 2 ) = 3 .
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The means, standard-deviations and the correlations are computed on all the observations where X are not missing (see Further Details).
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some pair of variables or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the eigenvectors and the eigenvalues.

On exit when LAST=false, FAILURE is always set to false.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, XEIGVEC contains the eigenvectors of the variances-covariances matrix, when LAST=true.
  • COV= false, XEIGVEC contains the eigenvectors of the correlation matrix, when LAST=true.

By default, the eigenvectors of the correlation matrix are output.

COV needs to be specified only on the last call to COMP_EOF_MISS (e.g., when LAST=true).

SORT (INPUT, OPTIONAL) character

Sort the eigenvalues into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The eigenvectors are reordered accordingly.

SORT needs to be specified only on the last call to COMP_EOF_MISS (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of an intermediate tridiagonal form T of the covariance matrix. The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(XEIGVAL). Convergence usually occurs in about 2 * size(XEIGVAL) QR sweeps.

The default is 30.

MAXITER needs to be specified only on the last call to COMP_EOF_MISS (e.g., when LAST=true).

XMEAN (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XMEAN is present, XMEAN contains the variable means computed from all non-missing observations in the data matrix.

XMEAN needs to be specified only on the last call to COMP_EOF_MISS (e.g., when LAST=true).

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

XSTD needs to be specified only on the last call to COMP_EOF_MISS (e.g., when LAST=true).

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

Further Details

The subroutine computes the Empirical Orthogonal Functions with only one pass through the data.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XEIGVAL, XEIGVEC, and XCORP(:,1) are globally set to XMISS.

The means and standard-deviations of the data matrix are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations. The eigenvectors and eigenvalues are computed from these bivariate statistics.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVAR) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

For more details on EOF or PCA analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_eof_miss2 ( x, first, last, xeigval, xcorp, xmiss, failure, dimvar, cov, maxiter, ortho, xmean, xstd, xeigvec )

Purpose

COMP_EOF_MISS2 computes estimates of Empirical Orthogonal Functions (EOF; also known as Principal Component Analysis) from a data matrix possibly containing missing values.

COMP_EOF_MISS2 computes all the eigenvalues, and optionally selected eigenvectors (by inverse iteration), of the covariance (or correlation matrix) from a data matrix.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Empirical Orthogonal Functions are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XEIGVAL (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_EOF_MISS2 (e.g., when FIRST=true), XEIGVAL contains temporary results from previous calls to COMP_EOF_MISS2. XEIGVAL should not be changed between calls to COMP_EOF_MISS2.

On exit, when LAST=true:

  • XEIGVAL(:,1) contains the eigenvalues of the variance-covariance (or correlation) matrix from the data matrix. The near zero eigenvalues made negative by round off errors or because the variance-covariance (or correlation) matrix from the data matrix with missing values is not positive definite are set to zero.
  • XEIGVAL(:,2) contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL(:,1).

The shape of XEIGVAL must verify:

  • size( XEIGVAL, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVAL, 2 ) = 2 .
XCORP (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_EOF_MISS2 (e.g., when FIRST=true), XCORP is used as workspace to accumulate quantities from previous calls to COMP_EOF_MISS2. XCORP should not be changed between calls to COMP_EOF_MISS2.

On exit:

  • XCORP(:,1) contains the correlation or variance-covariance matrix, as controlled by the COV argument, stored in symmetric storage mode. The upper triangle of the symmetric correlation or variance-covariance matrix is packed columnwise in the linear array XCORP(:,1). More precisely, the j-th column of this matrix is stored in the array XCORP(:,1) as follows:

    XCORP(i + (j-1) * j/2,1) = XCOR(i,j) for 1<=i<=j;

  • XCORP(:,2) contains the upper triangle of the matrix of the incidence values between each pair of variables, packed columnwise, in a linear array. XCORP(i + (j-1) * j/2,2) indicates the numbers of non-missing pairs which were used in the calculation of the covariance (or correlation) between variables i and j, for 1<=i<=j .

  • XCORP(:,3:4) is used as workspace and contains no useful informations.

The shape of XCORP must verify:

  • size( XCORP, 1 ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2 ,
  • size( XCORP, 2 ) = 4 .
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The means, standard-deviations and the correlations are computed on all the observations where X are not missing (see Further Details).
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some pair of variables or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the eigenvalues or that some eigenvectors failed to converge with MAXITER inverse iterations.

On exit when LAST=false, FAILURE is always set to false.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, when argument COV is present, COV is used as follows, if:

  • COV= true, the eigenvalues and the eigenvectors are computed from the variances-covariances matrix, when LAST=true.
  • COV= false, the eigenvalues and the eigenvectors are computed from the correlation matrix, when LAST=true.

By default, the eigenvalues and eigenvectors of the correlation matrix are computed.

COV needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine for computing the eigenvectors. By default, 2 inverse iterations are performed for all the eigenvectors. This optional argument is used only if the XEIGVEC is present.

MAXITER needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed eigenvectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the XEIGVEC is present.

The default is FALSE.

ORTHO needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

XMEAN (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XMEAN is present, XMEAN contains the variable means computed from all non-missing observations in the data matrix.

XMEAN needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XSTD is present, XSTD contains the variable standard-deviations.

XSTD needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XEIGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true, XEIGVEC contains the first size(XEIGVEC,2) eigenvectors of the variance-covariance (or correlation) matrix from the data matrix.

XEIGVEC needs to be specified only on the last call to COMP_EOF_MISS2 (e.g., when LAST=true).

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) <= size( X, DIMVAR ) .

Further Details

The subroutine computes the Empirical Orthogonal Functions with only one pass through the data.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XEIGVAL, XEIGVEC, and XCORP(:,1) are globally set to XMISS .

The means and standard-deviations of the data matrix are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations. The eigenvectors and eigenvalues are computed from these bivariate statistics.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVAR) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

For more details on EOF or PCA analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_eof_miss3 ( x, xmiss, dimvar, failure, xcorp, xincp, xeigval, xeigvec, maxiter, ortho )

Purpose

COMP_EOF_MISS3 computes estimates of Empirical Orthogonal Functions (EOF) or Principal Components (PC) from a data matrix X with n observations possibly containing missing values.

COMP_EOF_MISS3 computes an estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X, the associated matrix of incidence values, all the eigenvalues, and selected eigenvectors (by inverse iteration), of this estimate of the matrix product.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input data matrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables for which Empirical Orthogonal Functions or Principal Components are desired.
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) is computed on all the observations where X are not missing (see Further Details).
DIMVAR (INPUT) integer(i4b)

On entry, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input matrix X contains size(X,2) observations on size(X,1) variables and the matrix product (1/n) (X * X’) is computed.
  • DIMVAR = 2, the input matrix X contains size(X,1) observations on size(X,2) variables and the matrix product (1/n) (X’ * X) is computed.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than one valid observation were present for some pair of variables or that maximum accuracy was not achieved when computing the eigenvalues or that some eigenvectors failed to converge with MAXITER inverse iterations.
XCORP (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, XCORP(:) contains the estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’), stored in symmetric storage mode.

The upper triangle of the symmetric matrix product matrix is packed columnwise in the linear array XCORP(:). More precisely, the j-th column of this matrix is stored in the array XCORP(:) as follows:

XCORP(i + (j-1) * j/2,1) = XCOR(i,j) for 1<=i<=j;

The size of XCORP must verify: size( XCORP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2

XINCP (OUTPUT, OPTIONAL) integer(i4b), dimension(:)

On exit, XINCP(:) contains the upper triangle of the matrix of the incidence values between each pair of variables, packed columnwise, in a linear array. XINCP(i + (j-1) * j/2) indicates the numbers of non-missing pairs which were used in the calculation of the scalar product between variables i and j, for 1<=i<=j .

The size of XINCP must verify: size( XINCP ) = ( size(X,DIMVAR) * (size(X,DIMVAR)+1) )/2

XEIGVAL (OUTPUT, OPTIONAL) real(stnd), dimension(:,2)

On exit:

  • XEIGVAL(:,1) contains the eigenvalues in decreasing order of the estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X. The near zero eigenvalues made negative by round off errors or because the matrix product from the data matrix X with missing values is not positive definite are set to zero.
  • XEIGVAL(:,2) contains percentages of total variance associated with the eigenvectors in the order of the eigenvalues stored in XEIGVAL(:,1).

The shape of XEIGVAL must verify:

  • size( XEIGVAL, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVAL, 2 ) = 2 .
XEIGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, XEIGVEC contains the first size(XEIGVEC,2) eigenvectors of the estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) from the data matrix X.

The shape of XEIGVEC must verify:

  • size( XEIGVEC, 1 ) = size( X, DIMVAR ) ,
  • size( XEIGVEC, 2 ) <= size( X, DIMVAR ) .
MAXITER (INPUT, OPTIONAL) integer(i4b)
The number of inverse iterations performed in the subroutine for computing the eigenvectors. By default, 2 inverse iterations are performed for all the eigenvectors. This optional argument is used only if the XEIGVEC is present.
ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed eigenvectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the XEIGVEC is present.

The default is FALSE.

Further Details

The subroutine computes the Empirical Orthogonal Functions or the Principal Components with only one pass through the data.

The estimate of the matrix product (1/n) (X’ * X) or (1/n) (X * X’) is computed from all valid pairs of observations. The eigenvectors and eigenvalues are computed from these bivariate statistics.

If fewer than one valid observation is present for some pair of variables, the scalar product between this pair of variables is set to zero for computing the eigenvectors and eigenvalues of the estimated matrix product.

If size(X,3-DIMVAR)<=0, the subroutine set FAILURE to true and returns without doing any computations.

For more details on EOF or PCA analysis, see:

  1. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300

subroutine comp_pc_eof ( x, xeigvec, xsingval, xpc, dimvar, xmean, xstd, xpccor )

Purpose

COMP_PC_EOF computes estimates of Principal Components (PC) from a data matrix and a set of eigenvectors derived from an EOF or PCA analysis.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Principal Components are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
XEIGVEC (INPUT) real(stnd), dimension(:,:)

On entry, XEIGVEC contains selected eigenvectors of th variance-covariance (or correlation) matrix from the data matrix.

The shape of XEIGVEC must verify: size( XEIGVEC, 1 ) = size( X, DIMVAR ).

XSINGVAL (INPUT) real(stnd), dimension(:)

On entry, XSINGVAL must contain the singular values of the covariance (or correlation) matrix from the data matrix associated with the eigenvectors in XEIGVEC array. The Principal Components are normalized by XSINGVAL on output (the variances of the Principal Components are equal to one).

The size of XSINGVAL must verify: size( XSINGVAL ) = size( XEIGVEC, 2 ).

XPC (OUTPUT) real(stnd), dimension(:,:)

On exit, XPC contains the normalized Principal Components derived from X and XEIGVEC.

The shape of XPC must verify:

  • size( XPC, 1 ) = size( X, 3-DIMVAR ) ,
  • size( XPC, 2 ) = size( XEIGVEC, 2 ) .
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Principal Components are computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Principal Components are computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XPCCOR (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when XPCCOR is present, XPCCOR contains the correlations (or covariances) between the data matrix XX and the principal components (factor loadings matrix).

This optional argument may be specified in a call to COMP_PC_EOF with size(X,3-DIMVAR) = size(XPC,1) = 0 (e.g., with no observations) such as:

call comp_pc_eof( x(:nvar,1:0), xeigvec(:nvar,:neigvec), xsingval(:neigvec), &
xpc(1:0,:neigvec), xpccor=xpccor(:nvar,:neigvec) )

The shape of XPCCOR must verify:

  • size( XPCCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XPCCOR, 2 ) = size( XEIGVEC, 2 ) .

Further Details

The subroutine computes the Principal Components with only one pass through the data.

If unnormalized PCs are desired, use argument XSINGVAL with all values set to one, however in this case, do not use argument XPCCOR.

subroutine comp_ortho_rot_eof ( fac, rot_fac, orot, std_rot_fac, failure, knorm, maxiter, w, delta )

Purpose

COMP_ORTHO_ROT_EOF performs an orthogonal rotation of a (partial) EOF model (e.g., a factor loading matrix) using a generalized orthomax criterion, including quartimax, varimax and equamax rotation methods.

Arguments

FAC (INPUT) real(stnd), dimension(:,:)

On entry, the unrotated EOF model (e.g., the input factor loading matrix). The number of EOFs or factors in the model is equal to nf=size(FAC,2) and the number of variables is equal to nv=size(FAC,1).

The shape of FAC must verify:

  • size( FAC, 1 ) = nv >= size( FAC, 2 ) = nf >= 2 .
ROT_FAC (OUTPUT) real(stnd), dimension(:,:)

On exit, the rotated EOF model (e.g., the rotated factor loading matrix).

The shape of ROT_FAC must verify:

  • size( ROT_FAC, 1 ) = size( FAC, 1 ) = nv ,
  • size( ROT_FAC, 2 ) = size( FAC, 2 ) = nf .
OROT (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed orthogonal rotation matrix.

The shape of OROT must verify:

  • size( OROT, 1 ) = size( OROT, 2 ) = size( FAC, 2 ) = nf .
STD_ROT_FAC (OUTPUT) real(stnd), dimension(:)

On exit, the standard-deviations accounted for by the rotated EOFs.

The size of STD_ROT_FAC must verify:

  • size( STD_ROT_FAC ) = size( FAC, 2 ) = nf.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that convergence did not occur in MAXITER iterations. But, convergence was assumed and calculations continued so that the results can still be useful.
KNORM (INPUT, OPTIONAL) logical(lgl)

On entry, if:

  • KNORM = true indicates that the rows of the input unrotated EOF model must be normalized following Kaiser’s method.
  • KNORM = false indicates that row normalization is not required.

The default is KNORM=true.

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of iterations allowed for rotations. MAXITER <= 10 defaults to 10 iterations.

The default is MAXITER = 30.

W (INPUT, OPTIONAL) real(stnd)

Input constant factor used to define the method of rotation. W can be any positive real number, but best values lie between 0. and 5.*nf. W <= 0. defaults to 0. .

On input, if:

  • W = 0. the quartimax method is used.
  • W = 1. the varimax method is used.
  • W = size(FAC,2)/2. the equamax method is used.

See Further Details for more information.

The default is W = 1., which means that the varimax method of rotation is used by default.

DELTA (INPUT, OPTIONAL) real(stnd)

Input convergence constant for rotation (e.g., criterion function; see Further Details). When the relative change in the criterion function is less than DELTA from one iteration to the next, convergence is assumed. DELTA=0.001 is typical. DELTA <= epsilon(DELTA) defaults to epsilon(DELTA). Similarly, DELTA > 0.02 defaults to 0.02.

The default is DELTA = 0.001.

Further Details

The subroutine performs an orthogonal rotation according to a generalized orthomax criterion. In this analytic method of orthogonal rotation, a criterion function is defined as

Q = ( [sum i=1 to nv][sum j=1 to nf] ROT_FAC(i,j)**4 ) - (W/nv).( [sum j=1 to nf] ( [sum i=1 to nv] ROT_FAC(i,j)**2 )**2 )/nv

, where W is a positive user-specified constant yielding a family of rotations, nv is the number of variables and nf is the number of factors or EOFs and this function is maximized by finding a nf-by-nf orthogonal rotation matrix OROT such that

ROT_FAC = FAC * OROT

is a maximum of Q (here FAC is the specified matrix of the unrotated EOFs, e.g., the unrotated factor loading matrix).

W is a parameter determining the kind of solution to be computed and may be set as follows:

  • W = 0. is the quartimax method, which attempts to get each variable to load highly on only one (or a few) EOFs or factors.
  • W = 1. is the varimax method, which attempts to load highly a relatively low number of variables on each EOF or factor. Varimax is the most widely used method of orthogonal rotation.
  • W = nf/2. is the equamax method, which is a compromise of the above two.

W can be any positive real number, but best values lie in the closed interval [0., 5.*nf]. Generally, the larger W is, the more equal is the dispersion of the variance accounted for across the rotated factors.

The method for optimizing Q proceeds by accumulating simple rotations where a simple rotation is defined to be one in which Q is optimized for two columns of ROT_FAC and for which the requirement that ROT be an orthogonal matrix is satisfied. A single iteration is defined to be such that each of the nf.(nf-1) possible simple rotations is performed (where nf is the number of EOFs or factors).

When the relative change in the criterion function Q from one iteration to the next is less than DELTA (the user-specified convergence criterion), the algorithm stops. DELTA=0.001 is generaly sufficient. Alternatively, the algorithm stops when the user-specified maximum number of iterations, MAXITER, is reached. MAXITER=30 is usually sufficient.

Kaiser (row) normalization can be performed on the EOFs (e.g., the factor loadings) prior to the rotation via the optional logical paramter KNORM. If, on input KNORM=true, the rows of FAC are first “normalized” by dividing each row by the square root of the sum of its squared elements. After the rotation is complete, each row of ROT_FAC is “denormalized” by multiplication by its initial normalizing constant.

The documentation of this subroutine is partially adapted from the FROTA subroutine in the IMSL library, which performs exactly the same task.

For more details on orthogonal rotations for EOF or PCA analysis, see:

  1. Jackson, J.E., 2003:
    A user’s guide to principal components. John Wiley and Sons, New York, USA, Chapter 8, 592 pp., ISBN:978-0-471-47134-9
  2. Jolliffe, I.T., 2002:
    Principal component analysis. Springer-Verlag, New York, USA, Chapters 7 and 11, 487 pp., 2nd Ed, ISBN:978-0-387-22440-4
  3. von Storch, H., and Zwiers, F.W., 2002:
    Statistical Analysis in Climate Research. Cambridge, UK, Chapter 13, 484 pp., ISBN:9780521012300
  4. Jennrich, R.I., 1970:
    Orthogonal rotation algorithms. Psychometrika, Vol. 35, 229-235
  5. Clarkson, D.B., and Jennrich, R.I., 1988:
    Quartic rotation criteria and algorithms. Psychometrika, Vol. 53, 251-259
  6. Kaiser, H.F., 1958:
    The varimax criterion for analytic rotation in factor analysis. Psychometrika, Vol. 23, 187-200
  7. Jennrich, R.I., 2001:
    A simple general procedure for orthogonal rotation. Psychometrika, Vol. 66, 289-306
  8. Bernaards, C.A., and Jennrich, R.I., 2005:
    Gradient Projection Algorithms and Software for Arbitrary Rotation Criteria in Factor Analysis. Educational and Psychological Measurement, Vol. 65, 676-696

subroutine comp_smooth_rot_pc ( pc, std_pc, rot_pc, orot, std_rot_pc, failure, maxiter, d, smooth )

Purpose

COMP_SMOOTH_ROT_PC performs an orthogonal rotation of a (partial) EOF model (e.g., the standardized Principal Component time series) by minimizing a smoothness criterion.

Arguments

PC (INPUT) real(stnd), dimension(:,:)

On entry, the unrotated EOF model (e.g., the standardized Principal Component time series). The number of EOFs or factors in the model is equal to nf=size(PC,2) and the number of observations is equal to nobs=size(PC,1). The time observations in the original dataset must be ordered, but not neccessarily equally spaced. See Further details and the description the optional real vector argument D below for more information.

The shape of PC must verify:

  • size( PC, 1 ) = nobs >= size( PC, 2 ) = nf >= 2 .
STD_PC (INPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the input standardized Principal Component time series.

The size of STD_PC must verify:

  • size( STD_PC ) = size( PC, 2 ) = nf.
ROT_PC (OUTPUT) real(stnd), dimension(:,:)

On exit, the rotated EOF model (e.g., the rotated standardized Principal Component time series).

The shape of ROT_PC must verify:

  • size( ROT_PC, 1 ) = size( PC, 1 ) = nobs ,
  • size( ROT_PC, 2 ) = size( PC, 2 ) = nf .
OROT (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed orthogonal rotation matrix.

The shape of OROT must verify:

  • size( OROT, 1 ) = size( OROT, 2 ) = size( PC, 2 ) = nf .
STD_ROT_PC (OUTPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the rotated standardized Principal Component time series.

The size of STD_ROT_PC must verify:

  • size( STD_ROT_PC ) = size( PC, 2 ) = nf.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that convergence did not occur in MAXITER iterations. But, convergence was assumed and calculations continued so that the results can still be useful.
MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the symmetric matrix used to find the minima of the smoothness criterion.

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(PC,2). Convergence usually occurs in about 2 * size(PC,2) QR sweeps.

The default is 30.

D (INPUT, OPTIONAL) real(stnd), dimension(:)

Input vector indexing the time observations, if the time interval between sucessive observations is not constant. This optional argument has no effect on the computed solutions if the time interval is constant.

See Further Details for more information.

The size of D must verify the relation:

  • size( D ) = nobs .
SMOOTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

Smoothness criteria for the rotated standardized Principal Component time series.

See Further Details for more information.

The size of SMOOTH must verify the relation:

  • size( SMOOTH ) = nf .

Further Details

The subroutine performs an orthogonal rotation of standardized Principal Component time series according to a smoothness criterion, which is defined as

SM(Y) = || Y(i) - Y(i-1) ||**2 = [sum i=2 to nobs] (Y(i) - Y(i-1))**2

for a vector Y of nobs observations and satisfying || Y || = 1 (or Var(Y) = 1 ). More precisely, the subroutine finds the linear combinations of the input Principal Component time series, which give successively the minimum of SM(Y) with the requirement that the rotated Principal Component time series remain orthogonal (e.g., uncorrelated) and of norm unity. This is equivalent to find a nf-by-nf orthogonal rotation matrix OROT such that

ROT_PC = PC * OROT

and [sum i=1 to nf] SM( ROT_PC(:,i) ) is a minimum.

The smoothness criterion weights each squared successive difference equally, based on the assumption that the time observations are ordered and equally spaced. If the observations are not taken in successive, equally spaced time periods, the smoothness criterion can be suitably modified as

SM(Y) = [sum i=2 to nobs] ( (Y(i) - Y(i-1))/(D(i) - D(i-1)) )**2

where the D vector gives the spacing between observations, and so that each squared successive difference is weighted accordingly to the unequal spacing between observations when computing the smoothness criterion.

For more details on orthogonal rotations to smooth functions, see:

  1. Arbuckle, J., and Friendly, M.L., 1977:
    On rotating to smooth functions. Psychometrika, Vol. 42, 127-140
  2. Solow, A.R., and Patwardhan, A., 1996:
    Extracting a smooth trend from a time series: A modification of Singular Spectrum Analysis. Journal of Climate, Vol. 9, 2163-2166
  3. Jolliffe, I.T., 2002:
    Principal component analysis. Springer-Verlag, New York, USA, Chapters 7 and 11, 487 pp., 2nd Ed, ISBN:978-0-387-22440-4

subroutine comp_lfc_rot_pc ( pc, std_pc, nt, rot_pc, orot, std_rot_pc, failure, maxiter, itdeg, ntjump, residual, smooth )

Purpose

COMP_LFC_ROT_PC performs an orthogonal rotation of a (partial) EOF model (e.g., the standardized Principal Component time series) towards low-frequency or high-frequency components using the eigenvectors of the covariance matrix between the standardized Principal Component time series filtered with a LOESS smoother.

Arguments

PC (INPUT) real(stnd), dimension(:,:)

On entry, the unrotated EOF model (e.g., the standardized Principal Component time series). The number of EOFs or factors in the model is equal to nf=size(PC,2) and the number of observations is equal to nobs=size(PC,1). The time observations in the original dataset must be ordered and equally spaced.

The shape of PC must verify:

  • size( PC, 1 ) = nobs >= size( PC, 2 ) = nf >= 2 .
STD_PC (INPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the input standardized Principal Component time series.

The size of STD_PC must verify:

  • size( STD_PC ) = size( PC, 2 ) = nf.
NT (INPUT) integer(i4b)
On entry, the length of the LOESS trend smoother. The value of NT should be an odd integer greater than or equal to 3. As NT increases the values of the rotated standardized PC components become smoother if residual=false or rougher (e.g., high-frequency) if residual=true.
ROT_PC (OUTPUT) real(stnd), dimension(:,:)

On exit, the rotated EOF model (e.g., the rotated standardized Principal Component time series).

The shape of ROT_PC must verify:

  • size( ROT_PC, 1 ) = size( PC, 1 ) = nobs ,
  • size( ROT_PC, 2 ) = size( PC, 2 ) = nf .
OROT (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed orthogonal rotation matrix.

The shape of OROT must verify:

  • size( OROT, 1 ) = size( OROT, 2 ) = size( PC, 2 ) = nf .
STD_ROT_PC (OUTPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the rotated standardized Principal Component time series.

The size of STD_ROT_PC must verify:

  • size( STD_ROT_PC ) = size( PC, 2 ) = nf.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that convergence did not occur in MAXITER iterations. But, convergence was assumed and calculations continued so that the results can still be useful.
MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the covariance matrix used to find the orthogonal rotation matrix OROT.

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(PC,2). Convergence usually occurs in about 2 * size(PC,2) QR sweeps.

The default is 30.

ITDEG (INPUT, OPTIONAL) integer(i4b)

On entry, the degree of locally-fitted polynomial in LOESS trend smoothing. The value must be 0, 1 or 2.

By default, ITDEG is set to 1.

NTJUMP (INPUT, OPTIONAL) integer(i4b)

On entry, the skipping value for LOESS trend smoothing.

By default, NTJUMP is set to NT/10.

RESIDUAL (INPUT, OPTIONAL) logical(lgl)

On entry, if:

  • RESIDUAL = true : the rotation is done towards high-frequency components, e.g., using the eigenvectors of the covariance matrix between the residual from the trends of the original standardized Principal Component time series.
  • RESIDUAL = false : the rotation is done towards low-frequency components, e.g., using the eigenvectors of the covariance matrix between the trends of the original standardized Principal Component time series.

In both cases, the trends are estimated with a LOESS smoother determined by the values of the NT, ITDEG and NTJUMP arguments.

By default, RESIDUAL = false.

SMOOTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a vector containing the ratios of filtered variance to total variance for the rotated standardized Principal Component time series.

See Further Details for more information.

The size of SMOOTH must verify the relation:

  • size( SMOOTH ) = nf .

Further Details

The subroutine performs an orthogonal rotation of standardized Principal Component time series towards low- or high-frequency components using a framework described in the reference (1) as Low-Frequency Components Analysis (LFCA).

Here, LFCA is considered as an orthogonal rotation of the original Principal Component time series of a (partial) EOF model. The nf-by-nf orthogonal matrix used in the rotation is computed as the eigenvectors of the covariance (e.g., symmetric positive-definite) matrix between the filtered original standardized Principal Component time series. The filtering of the Principal Component time series is performed with the help of a LOESS smoother specified by the values of the NT, ITDEG and NTJUMP arguments.

Depending on the value of the optional logical argument RESIDUAL, the filtered standardized Principal Component time series are computed as residuals from the trends estimated with the LOESS smoother (if RESIDUAL = true) or as the trends estimated with the LOESS smoother (if RESIDUAL = false). In the first case, the rotated standardized Principal Component time series will be ordered from high-frequency to low-frequency modes and in the second case, they will be ordered from low-frequency to high-frequency modes.

For more details on orthogonal rotations to smooth functions, LFCA or LOESS smoothing, see:

  1. Wills, R.C., Schneider, T., Wallace, J.M., Battisti, D.S., and Hartmann, D.L., 2018:
    Disentangling Global Warming, Multidecadal Variability, and El Nino in Pacific Temperatures. Geophysical Research Letters, Vol. 45, 2487-2496
  2. Cleveland, W.S., 1979:
    Robust Locally Weighted Regression and Smoothing Scatterplots. Journal of the American Statistical Association, Vol. 74, 829-836
  3. Cleveland, W.S., and Devlin, S.J., 1988:
    Locally Weighted Regression: An Approach to Regression Analysis by Local Fitting. Journal of the American Statistical Association, Vol. 83, 596-610
  4. Cleveland, R.B., Cleveland, W.S., McRae, J.E., and Terpenning, I., 1990:
    STL: A Seasonal-Trend Decomposition Procedure Based on Loess. J. Official Stat., 6, 3-73

subroutine comp_filt_rot_pc ( pc, std_pc, pl, ph, rot_pc, orot, std_rot_pc, failure, maxiter, trend, win, smooth )

Purpose

COMP_FILT_ROT_PC performs an orthogonal rotation of a (partial) EOF model (e.g., the standardized Principal Component time series) towards low-frequency, high-frequency or band-pass components using the eigenvectors of the covariance matrix between the standardized Principal Component time series filtered with a windowed FFT filter.

Arguments

PC (INPUT) real(stnd), dimension(:,:)

On entry, the unrotated EOF model (e.g., the standardized Principal Component time series). The number of EOFs or factors in the model is equal to nf=size(PC,2) and the number of observations is equal to nobs=size(PC,1). The time observations in the original dataset must be ordered and equally spaced.

The shape of PC must verify:

  • size( PC, 1 ) = nobs >= size( PC, 2 ) = nf >= 2 .
STD_PC (INPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the input standardized Principal Component time series.

The size of STD_PC must verify:

  • size( STD_PC ) = size( PC, 2 ) = nf.
PL (INPUT) integer(i4b)

Minimum period of oscillation of desired components in number of timesteps for the filtered standardized PC time series used for computing the orthogonal rotation matrix.

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH,

PL must be equal to 0 or greater or equal to 2. Moreover, PL must be less or equal to nobs.

PH (INPUT) integer(i4b)

Maximum period of oscillation of desired components in number of timesteps for the filtered standardized PC time series used for computing the orthogonal rotation matrix.

USE PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

PH must be equal to 0 or greater or equal to 2. Moreover, PH must be less or equal to nobs.

ROT_PC (OUTPUT) real(stnd), dimension(:,:)

On exit, the rotated EOF model (e.g., the rotated standardized Principal Component time series).

The shape of ROT_PC must verify:

  • size( ROT_PC, 1 ) = size( PC, 1 ) = nobs ,
  • size( ROT_PC, 2 ) = size( PC, 2 ) = nf .
OROT (OUTPUT) real(stnd), dimension(:,:)

On exit, the computed orthogonal rotation matrix.

The shape of OROT must verify:

  • size( OROT, 1 ) = size( OROT, 2 ) = size( PC, 2 ) = nf .
STD_ROT_PC (OUTPUT) real(stnd), dimension(:)

On entry, the standard-deviations accounted for by the rotated standardized Principal Component time series.

The size of STD_ROT_PC must verify:

  • size( STD_ROT_PC ) = size( PC, 2 ) = nf.
FAILURE (OUTPUT) logical(lgl)

On exit:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that convergence did not occur in MAXITER iterations. But, convergence was assumed and calculations continued so that the results can still be useful.
MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the Schur decomposition of the covariance matrix used to find the orthogonal rotation matrix OROT.

The algorithm fails to converge if the number of QR sweeps exceeds MAXITER * size(PC,2). Convergence usually occurs in about 2 * size(PC,2) QR sweeps.

The default is 30.

TREND (INPUT, OPTIONAL) integer(i4b)

If:

  • TREND=+/-1 The means of the PC time series are removed before time filtering
  • TREND=+/-2 The drifts from the PC time series are removed before time filtering by using the formula: drift(:) = (PC(nobs,:nf) - PC(1,:nf))/(nobs - 1)
  • TREND=+/-3 The least-squares lines from the PC time series are removed before time filtering.

IF TREND=-1,-2 or -3, the means, drifts or least-squares lines are reintroduced post-filtering, respectively.

For other values of TREND nothing is done before or after filtering the standardized Principal Component time series.

WIN (INPUT, OPTIONAL) real(stnd)

By default, Hamming window filtering is used (i.e. WIN=0.54). SET WIN=0.5 for Hanning window or WIN=1 for rectangular window.

WIN must be greater or equal to O.5 and less or equal to 1.

SMOOTH (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, a vector containing the ratios of filtered variance to total variance for the rotated standardized Principal Component time series.

The size of SMOOTH must verify the relation:

  • size( SMOOTH ) = nf .

Further Details

The subroutine performs an orthogonal rotation of standardized Principal Component time series towards low-, high-frequency or band-pass components using a framework described in the reference (1) as Low-Frequency Components Analysis (LFCA).

Here, LFCA is considered as an orthogonal rotation of the original Principal Component time series of a (partial) EOF model. The nf-by-nf orthogonal matrix used in the rotation is computed as the eigenvectors of the covariance (e.g., symmetric positive-definite) matrix between the filtered original standardized Principal Component time series. The filtering of the Principal Component time series is performed with the help of a windowed FFT filter specified by the values of the PL, PH, TREND and WIN arguments. See reference (2) for more detailed on the windowed FFT filter used here. This windowed filter is also implemented in subroutine HWFILTER, which is included in module Time_Series_Procedures.

Use PL=0 for high-pass filtering frequencies corresponding to periods shorter than PH, or PH=0 for low-pass filtering frequencies corresponding to periods longer than PL.

Setting PH<PL is also allowed and performs band rejection of periods between PH and PL (i.e. in that case the meaning of the PL and PH arguments are reversed).

Examples:

  • For quarterly data, PL=6, PH=32 perform rotation towards Principal Component time series with periods between 1.5 and 8 yrs.
  • For monthly data, PL=0, PH=24 perform rotation towards Principal Component time series with periods less than 2 yrs.

Thus, depending on the values of the PL and PH arguments, the rotated standardized Principal Component time series will be ordered from high-frequency to low-frequency modes or vice-versa.

For more details on orthogonal rotations to smooth functions, LFCA or windowed filtering, see:

  1. Wills, R.C., Schneider, T., Wallace, J.M., Battisti, D.S., and Hartmann, D.L., 2018:
    Disentangling Global Warming, Multidecadal Variability, and El Nino in Pacific Temperatures. Geophysical Research Letters, Vol. 45, 2487-2496
  2. Iacobucci, A., and Noullez, A., 2005:
    A Frequency Selective Filter for Short-Length Time Series. Computational Economics, 25,75-102.

subroutine comp_mca ( x, y, first, last, xstat, ystat, xysingval, xsingvec, failure, dimvarx, dimvary, cov, sort, maxiter, ysingvec, xysingvar, xycor )

Purpose

COMP_MCA performs Maximum Covariance Analysis (MCA) or canonical covariance analysis between two data matrices XX and YY.

COMP_MCA computes the singular value decomposition (SVD) of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY. This SVD is written

XYCOR = U * SIGMA * V’

where SIGMA is a m-by-n matrix which is zero except for its min(m,n) diagonal elements, U is a m-by-m orthogonal matrix, and V is a n-by-n orthogonal matrix. The diagonal elements of SIGMA are the singular values of XYCOR; they are real and non-negative. The first min(m,n) columns of U and V are the left and right singular vectors of XYCOR.

The routine returns the singular values, the left and, optionally, the right singular vectors of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY.

Arguments

X (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(X,3-DIMVARX) observations on size(X,DIMVARX) variables from the “left” matrix of data XX. By default, DIMVARX is equal to 1. See description of optional DIMVARX argument for details. If all the data are available at once, X can be the full data matrix XX.

The shape of X must verify: size( X, 3-DIMVARX ) = size( Y, 3-DIMVARY ) .

Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the “right” matrix of data YY. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVARX ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA. XSTAT should not be changed between calls to COMP_MCA.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the “left” data matrix XX.
  • XSTAT(:,2) contains the standard-deviations of the “left” data matrix XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVARX ) ,
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA. YSTAT should not be changed between calls to COMP_MCA.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the “right” data matrix YY.
  • YSTAT(:,2) contains the the standard-deviations of the “right” data matrix YY.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 2 .
XYSINGVAL (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), XYSINGVAL(1) contains count of observations from previous calls to COMP_MCA. XYSINGVAL(1) should not be changed between calls to COMP_MCA.

On exit, XYSINGVAL contains the singular values of the correlation (or covariance) matrix XYCOR between the data matrices XX and YY.

The size of XYSINGVAL must verify: size( XYSINGVAL ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

XSINGVEC (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), XSINGVEC is used as workspace to accumulate quantities on previous calls to COMP_MCA. XSINGVEC should not be changed between calls to COMP_MCA.

On exit, when LAST=true, XSINGVEC is overwritten with the first min(size(X,DIMVARX),size(Y,DIMVARY)) columns of U, the left singular vectors of the correlation (or covariance) matrix XYCOV between XX and YY.

The shape of XSINGVEC must verify:

  • size( XSINGVEC, 1 ) = size( X, DIMVARX ) ,
  • size( XSINGVEC, 2 ) = size( Y, DIMVARY ).
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that maximum accuracy was not achieved when computing the SVD of the covariance (or correlation) matrix between the data matrices XX and YY .

On exit when LAST=false, FAILURE is always set to false.

DIMVARX (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARX is present, DIMVARX is used as follows, if:

  • DIMVARX = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVARX = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVARX = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows, if:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV is present and COV=true, a covariance matrix between the data matrices XX and YY is computed instead of a correlation matrix.

By default, a correlation matrix is computed.

COV needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

SORT (INPUT, OPTIONAL) character

Sort the singular values into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The singular vectors are rearranged accordingly.

SORT needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of XYCOR fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default is 10.

MAXITER needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

YSINGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true and YSINGVEC is present, YSINGVEC contains the first min(size(X,DIMVARX),size(Y,DIMVARY)) columns of V, the right singular vectors of the correlation (or covariance) matrix XYCOR between XX and YY.

YSINGVEC needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

The shape of YSINGVEC must verify:

  • size( YSINGVEC, 1 ) = size( Y, DIMVARY ) ,
  • size( YSINGVEC, 2 ) = min( size(X,DIMVARX), size(Y,DIMVARY) ) .
XYSINGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XYSINGVAR is present, XYSINGVAR contains percentages of total squared covariance associated with the left and right singular vectors in order of the singular values stored in XYSINGVAL.

XYSINGVAR needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

The size of XYSINGVAR must verify: size( XYSINGVAR ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

XYCOR (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true and XYCOR is present, XYCOR contains the correlation or variance-covariance matrix between data arrays XX and YY, as controlled by the COV argument.

XYCOR needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVARX ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) .

Further Details

The subroutine computes the basic univariate statistics and the correlation (or covariance) matrix with only one pass through the data.

If fewer than two valid observations were present, the statistics XSTAT, YSTAT, XYSINGVAL XSINGVEC, YSINGVEC, XYSINGVAR and XYCOR are set to Nan code.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVARX) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

subroutine comp_mca2 ( x, y, first, last, xstat, ystat, xysingval, xycor, failure, dimvarx, dimvary, cov, savecor, maxiter, ortho, xysingvar, xysingvec )

Purpose

COMP_MCA2 performs Maximum Covariance Analysis (MCA) or canonical covariance analysis between two data matrices XX and YY.

COMP_MCA2 computes a partial singular value decomposition (SVD) of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY. This partial SVD is written

U(:m,:k) * SIGMA(:k,:k) * V(:n,:k)’

where SIGMA is a k-by-k matrix which is zero except for its k diagonal elements, U is a m-by-k orthogonal matrix, and V is a n-by-k orthogonal matrix. The diagonal elements of SIGMA are the first k singular values of XYCOR; they are real and non-negative. The k columns of U and V are the first k left and right singular vectors of XYCOR.

COMP_MCA2 computes all the singular values, and, optionally, selected left and right singular vectors (by inverse iteration), of the covariance (or correlation matrix) XYCOR between two data matrices XX and YY.

Arguments

X (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(X,3-DIMVARX) observations on size(X,DIMVARX) variables from the “left” matrix of data XX. By default, DIMVARX is equal to 1. See description of optional DIMVARX argument for details. If all the data are available at once, X can be the full data matrix XX.

The shape of X must verify: size( X, 3-DIMVARX ) = size( Y, 3-DIMVARY ) .

Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the “right” matrix of data YY. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVARX ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA. XSTAT should not be changed between calls to COMP_MCA.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the “left” data matrix XX.
  • XSTAT(:,2) contains the standard-deviations of the “left” data matrix XX.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVARX ) ,
  • size( XSTAT, 2 ) = 2 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,2)

On entry, after the first call to COMP_MCA (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA. YSTAT should not be changed between calls to COMP_MCA.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the “right” data matrix YY.
  • YSTAT(:,2) contains the the standard-deviations of the “right” data matrix YY.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 2 .
XYSINGVAL (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_MCA2 (e.g., when FIRST=true), XYSINGVAL(1) contains count of observations from previous calls to COMP_MCA2. XYSINGVAL(1) should not be changed between calls to COMP_MCA2.

On exit, XYSINGVAL contains the singular values of the correlation (or covariance) matrix XYCOR between the data matrices XX and YY.

The size of XYSINGVAL must verify: size( XYSINGVAL ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_MCA2 (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_MCA2. XYCOR should not be changed between calls to COMP_MCA2.

On exit, when LAST=true and SAVECOR=true, XYCOR contains the correlation or variance-covariance matrix as controlled by the COV argument. In this case XYCOR(i,j) contains the correlation (or covariance) coefficient between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVARX=2 and DIMVARY=2 ).

If SAVECOR=false, the correlation (or covariance) matrix is not saved on exit. In this case, XYCOR does not contain useful information.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVARX ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) .
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present or that maximum accuracy was not achieved when computing the singular values or that some singular vectors failed to converge with MAXITER inverse iterations.

On exit when LAST=false, FAILURE is always set to false.

DIMVARX (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARX is present, DIMVARX is used as follows, if:

  • DIMVARX = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVARX = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVARX = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows, if:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV is present and COV=true, a covariance matrix between the data matrices XX and YY is computed instead of a correlation matrix.

By default, a correlation matrix is computed.

COV needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

SAVECOR (INPUT, OPTIONAL) logical(lgl)

On exit, when argument SAVECOR is present and LAST=true, SAVECOR is used as follows, if:

  • SAVECOR= true, the correlation (or covariance) matrix is saved in argument XYCOR.
  • SAVECOR= false, the correlation (or covariance) matrix is destroyed.

By default, the correlation (or covariance) matrix is destroyed.

SAVECOR needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine for computing the singular vectors. By default, 2 inverse iterations are performed for all the singular vectors. This optional argument is used only if the XYSINGVEC argument is present.

MAXITER needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed singular vectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the XYSINGVEC argument is present.

The default is FALSE.

ORTHO needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

XYSINGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XYSINGVAR is present, XYSINGVAR contains percentages of total squared covariance associated with the left and right singular vectors in order of the singular values stored in XYSINGVAL.

XYSINGVAR needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

The size of XYSINGVAR must verify: size( XYSINGVAR ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

XYSINGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true and XYSINGVEC is present, XYSINGVEC contains the first columns of U and V, the first k left and right singular vectors of the correlation (or covariance) matrix XYCOR between XX and YY.

The first k left singular vectors are stored in XYSINGVEC(:size(X,DIMVARX),:) The first k right singular vectors are stored in XYSINGVEC(size(X,DIMVARX)+1:,:)

XYSINGVEC needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

The shape of XYSINGVEC must verify:

  • size( XYSINGVEC, 1 ) = size( X, DIMVARX ) + size( Y, DIMVARY ) ,
  • size( XYSINGVEC, 2 ) <= min( size(X,DIMVARX) ,size(Y,DIMVARY) ).

Further Details

The subroutine computes the basic univariate statistics and the correlation (or covariance) matrix with only one pass through the data.

If fewer than two valid observations were present, the statistics XSTAT, YSTAT, XYSINGVAL XYCOR , XYSINGVAR and XYSINGVEC are set to Nan code.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVARX) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

subroutine comp_mca_miss ( x, y, first, last, xstat, ystat, xycor, xymiss, failure, dimvarx, dimvary, cov, sort, maxiter, xysingval, xysingvar, ysingvec )

Purpose

COMP_MCA_MISS performs Maximum Covariance Analysis (MCA) or canonical covariance analysis between two data matrices XX and YY possibly containing missing values.

COMP_MCA_MISS computes the singular value decomposition (SVD) of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY. This SVD is written

XYCOR = U * SIGMA * V’

where SIGMA is a m-by-n matrix which is zero except for its min(m,n) diagonal elements, U is a m-by-m orthogonal matrix, and V is a n-by-n orthogonal matrix. The diagonal elements of SIGMA are the singular values of XYCOR; they are real and non-negative. The first min(m,n) columns of U and V are the left and right singular vectors of XYCOR.

The routine returns the singular values, the left and, optionally, the right singular vectors of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY.

Arguments

X (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(X,3-DIMVARX) observations on size(X,DIMVARX) variables from the “left” matrix of data XX. By default, DIMVARX is equal to 1. See description of optional DIMVARX argument for details. If all the data are available at once, X can be the full data matrix XX.

The shape of X must verify: size( X, 3-DIMVARX ) = size( Y, 3-DIMVARY ) .

Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the “right” matrix of data YY. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVARX ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_MCA_MISS (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS. XSTAT should not be changed between calls to COMP_MCA_MISS.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the “left” data matrix XX.
  • XSTAT(:,2) contains the standard-deviations of the “left” data matrix XX.
  • XSTAT(:,3) contains the the numbers of non-missing observations in the “left” data matrix XX.
  • XSTAT(:,4) is used as workspace.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVARX ) ,
  • size( XSTAT, 2 ) = 4 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_MCA_MISS (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS. YSTAT should not be changed between calls to COMP_MCA_MISS.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the “right” data matrix YY.
  • YSTAT(:,2) contains the the standard-deviations of the “right” data matrix YY.
  • YSTAT(:,3) contains the the numbers of non-missing observations in the “right” data matrix YY.
  • YSTAT(:,4) is used as workspace.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 4 .
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:,4)

On entry, after the first call to COMP_MCA_MISS (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS. XYCOR should not be changed between calls to COMP_MCA_MISS.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(i,j,1) contains the correlation coefficients between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVARX=2 and DIMVARY=2 ).
  • XYCOR(i,j,2) contains the incidence values between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVARX=2 and DIMVARY=2). XYCOR(i,j,2) indicates the numbers of non-missing pairs of observations which were used in the calculation of XYCOR(i,j,1).
  • XYCOR(:,:,3) contains the first min(size(X,DIMVARX),size(Y,DIMVARY)) columns of U, the left singular vectors of the correlation (or covariance) matrix XYCOV between XX and YY.
  • XYCOR(:,:,4) is used as workspace.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVARX ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) ,
  • size( XYCOR, 3 ) = 4 .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate statistics and correlation (or covariance) matrix are computed on all the observations where X and Y are not missing (see Further Details).
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some pair of variables or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the SVD of the covariance (or correlation) matrix between the data matrices XX and YY .

On exit when LAST=false, FAILURE is always set to false.

DIMVARX (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARX is present, DIMVARX is used as follows, if:

  • DIMVARX = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVARX = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVARX = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows, if:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV is present and COV=true, a covariance matrix between the data matrices XX and YY is computed instead of a correlation matrix.

By default, a correlation matrix is computed.

COV needs to be specified only on the last call to COMP_MCA_MISS (e.g., when LAST=true).

SORT (INPUT, OPTIONAL) character

Sort the singular values into ascending order if SORT = ‘A’ or ‘a’, or in descending order if SORT = ‘D’ or ‘d’. The singular vectors are rearranged accordingly.

SORT needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

MAXITER (INPUT, OPTIONAL) integer(i4b)

MAXITER controls the maximum number of QR sweeps in the bidiagonal SVD phase of the SVD algorithm. The bidiagonal SVD algorithm of an intermediate bidiagonal form B of XYCOR fails to converge if the number of QR sweeps exceeds MAXITER * min(m,n). Convergence usually occurs in about 2 * min(m,n) QR sweeps.

The default is 10.

MAXITER needs to be specified only on the last call to COMP_MCA (e.g., when LAST=true).

XYSINGVAL (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, XYSINGVAL contains the singular values of the correlation (or covariance) matrix XYCOR between the data matrices XX and YY. If this optional argument is absent, the singular values are stored in XSTAT(:,4), when LAST=true.

The size of XYSINGVAL must verify: size( XYSINGVAL ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

YSINGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true and YSINGVEC is present, YSINGVEC contains the first min(size(X,DIMVARX),size(Y,DIMVARY)) columns of V, the right singular vectors of the correlation (or covariance) matrix XYCOR between XX and YY.

YSINGVEC needs to be specified only on the last call to COMP_MCA_MISS (e.g., when LAST=true).

The shape of YSINGVEC must verify:

  • size( YSINGVEC, 1 ) = size( Y, DIMVARY ) ,
  • size( YSINGVEC, 2 ) = min( size(X,DIMVARX), size(Y,DIMVARY) ) .
XYSINGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XYSINGVAR is present, XYSINGVAR contains percentages of total squared covariance associated with the left and right singular vectors in order of the singular values stored in XYSINGVAL.

XYSINGVAR needs to be specified only on the last call to COMP_MCA_MISS (e.g., when LAST=true).

The size of XYSINGVAR must verify: size( XYSINGVAR ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

Further Details

The subroutine computes the basic univariate statistics and the correlation (or covariance) matrix with only one pass through the data.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XYSINGVAL, XYSINGVEC, XYSINGVAR and XYCOR(:,:,1) are globally set to XMISS.

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations. The singular vectors and singular values are computed from these bivariate statistics.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVARX) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

subroutine comp_mca_miss2 ( x, y, first, last, xstat, ystat, xycor, xymiss, failure, dimvarx, dimvary, cov, xysingval, maxiter, ortho, xysingvar, xysingvec )

Purpose

COMP_MCA_MISS2 performs Maximum Covariance Analysis (MCA) or canonical covariance analysis between two data matrices XX and YY possibly containing missing values.

COMP_MCA_MISS2 computes a partial singular value decomposition (SVD) of the correlation (or covariance) matrix XYCOR between two data matrices XX and YY. This partial SVD is written

U(:m,:k) * SIGMA(:k,:k) * V(:n,:k)’

where SIGMA is a k-by-k matrix which is zero except for its k diagonal elements, U is a m-by-k orthogonal matrix, and V is a n-by-k orthogonal matrix. The diagonal elements of SIGMA are the first k singular values of XYCOR; they are real and non-negative. The k columns of U and V are the first k left and right singular vectors of XYCOR.

COMP_MCA_MISS2 computes all the singular values, and, optionally, selected left and right singular vectors (by inverse iteration), of the covariance (or correlation matrix) XYCOR between two data matrices XX and YY.

Arguments

X (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(X,3-DIMVARX) observations on size(X,DIMVARX) variables from the “left” matrix of data XX. By default, DIMVARX is equal to 1. See description of optional DIMVARX argument for details. If all the data are available at once, X can be the full data matrix XX.

The shape of X must verify: size( X, 3-DIMVARX ) = size( Y, 3-DIMVARY ) .

Y (INPUT) real(stnd), dimension(:,:)

On entry, input submatrix containing size(Y,3-DIMVARY) observations on size(Y,DIMVARY) variables from the “right” matrix of data YY. By default, DIMVARY is equal to 1. See description of optional DIMVARY argument for details. If all the data are available at once, Y can be the full data matrix YY.

The shape of Y must verify: size( Y, 3-DIMVARY ) = size( X, 3-DIMVARX ) .

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrices are the first submatrices of the data matrices XX and YY.
  • FIRST = false the current submatrices are not the first submatrices of the data matrices XX and YY.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrices are the last submatrices of the data matrices XX and YY.
  • LAST = false the current submatrix are not the last submatrices of the data matrices XX and YY.
XSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_MCA_MISS2 (e.g., when FIRST=true), XSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS2. XSTAT should not be changed between calls to COMP_MCA_MISS2.

On exit, when LAST=true, XSTAT contains the following statistics on all variables from the XX matrix:

  • XSTAT(:,1) contains the mean values of the “left” data matrix XX.
  • XSTAT(:,2) contains the standard-deviations of the “left” data matrix XX.
  • XSTAT(:,3) contains the the numbers of non-missing observations in the “left” data matrix XX.
  • XSTAT(:,4) is used as workspace.

The shape of XSTAT must verify:

  • size( XSTAT, 1 ) = size( X, DIMVARX ) ,
  • size( XSTAT, 2 ) = 4 .
YSTAT (INPUT/OUTPUT) real(stnd), dimension(:,4)

On entry, after the first call to COMP_MCA_MISS2 (e.g., when FIRST=true), YSTAT is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS2. YSTAT should not be changed between calls to COMP_MCA_MISS2.

On exit, when LAST=true, YSTAT contains the following statistics on all variables from the YY matrix:

  • YSTAT(:,1) contains the mean values of the “right” data matrix YY.
  • YSTAT(:,2) contains the the standard-deviations of the “right” data matrix YY.
  • YSTAT(:,3) contains the the numbers of non-missing observations in the “right” data matrix YY.
  • YSTAT(:,4) is used as workspace.

The shape of YSTAT must verify:

  • size( YSTAT, 1 ) = size( Y, DIMVARY ) ,
  • size( YSTAT, 2 ) = 4 .
XYCOR (INPUT/OUTPUT) real(stnd), dimension(:,:,4)

On entry, after the first call to COMP_MCA_MISS2 (e.g., when FIRST=true), XYCOR is used as workspace to accumulate quantities on previous calls to COMP_MCA_MISS2. XYCOR should not be changed between calls to COMP_MCA_MISS2.

On exit, when LAST=true, XYCOR contains the following statistics:

  • XYCOR(i,j,1) contains the correlation coefficients between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVARX=2 and DIMVARY=2 ).
  • XYCOR(i,j,2) contains the incidence values between XX(i,:) and YY(j,:) ( XX(:,i) and YY(:,j) if DIMVARX=2 and DIMVARY=2). XYCOR(i,j,2) indicates the numbers of non-missing pairs of observations which were used in the calculation of XYCOR(i,j,1).
  • XYCOR(:,:,3:4) is used as workspace.

The shape of XYCOR must verify:

  • size( XYCOR, 1 ) = size( X, DIMVARX ) ,
  • size( XYCOR, 2 ) = size( Y, DIMVARY ) ,
  • size( XYCOR, 3 ) = 4 .
XYMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X or Y which is equal to XYMISS is assumed to be missing or invalid. The basic univariate statistics and correlation (or covariance) matrix are computed on all the observations where X and Y are not missing (see Further Details).
FAILURE (OUTPUT) logical(lgl)

On exit when LAST=true:

  • FAILURE = false: indicates successful exit.
  • FAILURE = true: indicates that fewer than two valid observations were present for some pair of variables or that the observations on some variable were constant and the correlations were requested or that maximum accuracy was not achieved when computing the eigenvalues or that some eigenvectors failed to converge with MAXITER inverse iterations.

On exit when LAST=false, FAILURE is always set to false.

DIMVARX (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARX is present, DIMVARX is used as follows, if:

  • DIMVARX = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVARX = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables, respectively.

The default is DIMVARX = 1.

DIMVARY (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVARY is present, DIMVARY is used as follows, if:

  • DIMVARY = 1, the input submatrix Y contains size(Y,2) observations on size(Y,1) variables.
  • DIMVARY = 2, the input submatrix Y contains size(Y,1) observations on size(Y,2) variables, respectively.

The default is DIMVARY = 1.

COV (INPUT, OPTIONAL) logical(lgl)

On entry, if COV is present and COV=true, a covariance matrix between the data matrices XX and YY is computed instead of a correlation matrix.

By default, a correlation matrix is computed.

COV needs to be specified only on the last call to COMP_MCA_MISS2 (e.g., when LAST=true).

XYSINGVAL (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, XYSINGVAL contains the singular values of the correlation (or covariance) matrix XYCOR between the data matrices XX and YY. If this optional argument is absent, the singular values are stored in XSTAT(:,4), when LAST=true.

The size of XYSINGVAL must verify: size( XYSINGVAL ) = min( size(X,DIMVARX), size(Y,DIMVARY) ) .

MAXITER (INPUT, OPTIONAL) integer(i4b)

The number of inverse iterations performed in the subroutine for computing the singular vectors. By default, 2 inverse iterations are performed for all the singular vectors. This optional argument is used only if the XYSINGVEC argument is present.

MAXITER needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

ORTHO (INPUT, OPTIONAL) logical(lgl)

If ORTHO=true the computed singular vectors are orthogonalized by the Modified Gram-Schmidt algorithm. This optional argument is used only if the XYSINGVEC argument is present.

The default is FALSE.

ORTHO needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

XYSINGVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XYSINGVAR is present, XYSINGVAR contains percentages of total squared covariance associated with the left and right singular vectors in order of the singular values stored in XYSINGVAL.

XYSINGVAR needs to be specified only on the last call to COMP_MCA_MISS2 (e.g., when LAST=true).

The size of XYSINGVAR must verify: size( XYSINGVAR ) = min( size(X,DIMVARX), size(Y,DIMVARY) ).

XYSINGVEC (OUTPUT, OPTIONAL) real(stnd), dimension(:,:)

On exit, when LAST=true and XYSINGVEC is present, XYSINGVEC contains the first columns of U and V, the first k left and right singular vectors of the correlation (or covariance) matrix XYCOR between XX and YY.

The first k left singular vectors are stored in XYSINGVEC(:size(X,DIMVARX),:) The first k right singular vectors are stored in XYSINGVEC(size(X,DIMVARX)+1:,:)

XYSINGVEC needs to be specified only on the last call to COMP_MCA2 (e.g., when LAST=true).

The shape of XYSINGVEC must verify:

  • size( XYSINGVEC, 1 ) = size( X, DIMVARX ) + size( Y, DIMVARY ) ,
  • size( XYSINGVEC, 2 ) <= min( size(X,DIMVARX) ,size(Y,DIMVARY) ) .

Further Details

The subroutine computes the basic univariate statistics and the correlation (or covariance) matrix with only one pass through the data.

If fewer than two valid observations were present for some pair of variables or if the observations on some variable were constants, the statistics XYSINGVAL, XYSINGVEC, XYSINGVAR and XYCOR(:,:,1) are globally set to XMISS.

The means and standard-deviations of XX and YY are computed from all valid data. The correlation coefficients are based on these univariate statistics and on all valid pairs of observations. The singular vectors and singular values are computed from these bivariate statistics.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVARX) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

subroutine comp_pc_mca ( x, xsingvec, first, last, xpccor, pccorp, xpc, xn, dimvar, xmean, xstd, xpcvar )

Purpose

COMP_PC_MCA computes estimates of Singular Variables (SV) and correlation (or covariance) fields from a data matrix XX and a set of singular vectors derived from MCA analysis of XX with another matrix YY.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data for which Singular Variables are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix.
XSINGVEC (INPUT) real(stnd), dimension(:,:)

On entry, XSINGVEC contains selected left (or right) singular vectors of the SVD of the covariance (or correlation) matrix between the data matrix XX and another data matrix YY.

The shape of XSINGVEC must verify: size( XSINGVEC, 1 ) = size( X, DIMVAR ).

FIRST (INPUT) logical(lgl)

On entry, if:

  • FIRST = true the current submatrix is the first submatrix of the data matrix.
  • FIRST = false the current submatrix is not the first submatrix of the data matrix.
LAST (INPUT) logical(lgl)

On entry, if:

  • LAST = true the current submatrix is the last submatrix of the data matrix.
  • LAST = false the current submatrix is not the last submatrix of the data matrix.
XPCCOR (INPUT/OUTPUT) real(stnd), dimension(:,:)

On entry, after the first call to COMP_PC_MCA (e.g., when FIRST=true), XPCCOR is used as workspace to accumulate quantities on previous calls to COMP_PC_MCA. XPCCOR should not be changed between calls to COMP_PC_MCA.

On exit, when LAST=true, XPCCOR contains:

  • the correlations between the data matrix and the Singular Variables if the optional arguments XMEAN and XSTD are present.
  • the covariances between the data matrix and the normalized Singular Variables if only the optional argument XMEAN is present.

The shape of XPCCOR must verify:

  • size( XPCCOR, 1 ) = size( X, DIMVAR ) ,
  • size( XPCCOR, 2 ) = size( XSINGVEC, 2 ) .
PCCORP (INPUT/OUTPUT) real(stnd), dimension(:)

On entry, after the first call to COMP_PC_MCA (e.g., when FIRST=true), PCCORP is used as workspace to accumulate quantities on previous calls to COMP_PC_MCA. PCCORP should not be changed between calls to COMP_PC_MCA.

On exit, when LAST=true, PCCORP contains the correlation matrix between the Singular Variables stored in argument XPC. PCCORP is stored in symmetric storage mode (see further details).

The size of PCCORP must verify: size( PCCORP ) = ( size(XSINGVEC,2) * (size(XSINGVEC,2)+1) )/2.

XPC (OUTPUT) real(stnd), dimension(:,:)

On exit, XPC contains the unnormalized Singular Variables derived from X and XSINGVEC.

The shape of XPC must verify:

  • size( XPC, 1 ) = size( X, 3-DIMVAR ) ,
  • size( XPC, 2 ) = size( XSINGVEC, 2 ) .
XN (INPUT/OUTPUT) real(stnd)

On entry, after the first call to COMP_PC_MCA (e.g., when FIRST=true), XN contains count of observations from previous calls to COMP_PC_MCA. XN should not be changed between calls to COMP_PC_MCA.

On exit, XN contains the number of observations in the data matrix XX.

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Singular Variables are computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Singular Variables are computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XPCVAR (OUTPUT, OPTIONAL) real(stnd), dimension(:)

On exit, when LAST=true and XPCVAR is present, XPCVAR contains the variances of the Singular Variables stored in argument XPC.

The size of XPCVAR must verify: size( XPCVAR ) = size( XSINGVEC, 2 )

XPCVAR needs to be specified only on the last call to COMP_PC_MCA (e.g., when LAST=true).

Further Details

The subroutine computes the Singular Variables and the correlations matrices with only one pass through the data.

On exit, the upper triangle of the symmetric correlation matrix COR between the Singular Variables is packed columnwise in the linear array PCCORP. More precisely, the j-th column of this matrix COR is stored in the array PCCORP as follows:

PCCORP(i + (j-1) * j/2) = COR(i,j) for 1<=i<=j;

If fewer than two valid observations were present, the arguments XPCVAR, XPCCOR and PCCORP are set to nan() code.

This subroutine may be used in a call with no observations (e.g., size(X,3-DIMVAR) = size(XPC,1) = 0) in order to finish the computations with LAST=true when the total number of observations is unknown at the beginning of the computations.

The correlations (or covariance) between the Singular variables XPC and the data array YY may be computed easily from the outputs of COMP_MCA and COMP_PC_MCA. If YSINGVEC(:p,:k) are the k singular vectors (derived from data array YY) associated with the k singular values XYSINGVAL(:k) and the k singular vectors XSINGVEC(:m,:k) (derived from data array XX). Then, the correlations (or covariance) between the Singular variables XPC and the data array YY are equal to

spread( XYSINGVAL(:k)/SQRT(XPCVAR(:k)), dim=1, ncopies=p) * YSINGVEC(:p,:k)

This matrix is a covariance matrix between data array YY and the normalized Singular Variables XPC if a covariance matrix between data arrays XX and YY is analysed or a correlation matrix between data array YY and the Singular Variables XPC if a correlation matrix between data arrays XX and YY is analysed.

subroutine comp_pc ( x, xeigvec, xpc, dimvar, xmean, xstd, xsingval )

Purpose

COMP_PC computes estimates of a Principal Component (PC) from a data matrix XX and an eigenvector or singular vector derived from EOF or MCA analysis.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which the Principal Component is desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
XEIGVEC (INPUT) real(stnd), dimension(:)

On entry, XEIGVEC contains an eigenvector of the variance-covariance (or correlation) matrix from the data matrix XX or a selected left (or right) singular vector of the SVD of the covariance matrix between the data matrix XX and another data matrix YY.

The shape of XEIGVEC must verify: size( XEIGVEC ) = size( X, DIMVAR ).

XPC (OUTPUT) real(stnd), dimension(:)

On exit, XPC contains the Principal Component derived from X and XEIGVEC.

The size of XPC must verify: size( XPC ) = size( X, 3-DIMVAR ).

DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Principal Component is computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Principal Component is computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSINGVAL (INPUT, OPTIONAL) real(stnd)
On entry, XSINGVAL must contain the singular value of the covariance (or correlation) matrix from the data matrix XX associated with the eigenvector in XEIGVEC array. If SINGVAL is present, the Principal Component is normalized by XSINGVAL on output (the variance of the Principal Component is equal to one). This optional argument is useful only if XEIGVEC contains an eigenvector derived from an EOF analysis.

Further Details

The subroutine computes the Principal Component with only one pass through the data.

subroutine comp_pc ( x, xeigvec, xpc, dimvar, xmean, xstd, xsingval )

Purpose

COMP_PC computes estimates of Principal Components (PC) from a data matrix XX and a set of eigenvectors or singular vectors derived from EOF or MCA analysis.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which Principal Components are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
XEIGVEC (INPUT) real(stnd), dimension(:,:)

On entry, XEIGVEC contains selected eigenvectors of the variance-covariance (or correlation) matrix from the data matrix XX or selected left (or right) singular vectors of the SVD of the covariance matrix between the data matrix XX and another data matrix YY.

The shape of XEIGVEC must verify: size( XEIGVEC, 1 ) = size( X, DIMVAR ).

XPC (OUTPUT) real(stnd), dimension(:,:)

On exit, XPC contains the Principal Components derived from X and XEIGVEC.

The shape of XPC must verify:

  • size( XPC, 1 ) = size( X, 3-DIMVAR ) ,
  • size( XPC, 2 ) = size( XEIGVEC, 2 ) .
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Principal Components are computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Principal Components are computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSINGVAL (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, XSINGVAL must contain the singular values of the covariance (or correlation) matrix from the data matrix XX associated with the eigenvectors in XEIGVEC array. If SINGVAL is present, the Principal Components are normalized by XSINGVAL on output (the variances of the Principal Components are equal to one). This optional argument is useful only if XEIGVEC contains eigenvectors derived from an EOF analysis.

The size of XSINGVAL must verify: size( XSINGVAL ) = size( XEIGVEC, 2 ).

Further Details

The subroutine computes the Principal Components with only one pass through the data.

subroutine comp_pc_miss ( x, xeigvec, xpc, xmiss, dimvar, xmean, xstd, xsingval )

Purpose

COMP_PC_MISS computes estimates of a Principal Component (PC) from a data matrix XX and an eigenvector or singular vector derived from EOF or MCA analysis when XX contains missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which a Principal Components is desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
XEIGVEC (INPUT) real(stnd), dimension(:)

On entry, XEIGVEC contains a selected eigenvector of the variance-covariance (or correlation) matrix from the data matrix XX or a selected left (or right) singular vector of the SVD of the covariance matrix between the data matrix XX and another data matrix YY.

The shape of XEIGVEC must verify: size( XEIGVEC ) = size( X, DIMVAR ).

XPC (OUTPUT) real(stnd), dimension(:)

On exit, XPC contains the Principal Component derived from X and XEIGVEC.

The shape of XPC must verify: size( XPC ) = size( X, 3-DIMVAR ).

XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The Principal Components are computed from all the observations where X X are not missing by regression (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Principal Component is computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Principal Component is computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSINGVAL (INPUT, OPTIONAL) real(stnd)
On entry, XSINGVAL must contain the singular value of the covariance (or correlation) matrix from the data matrix XX associated with the eigenvector in XEIGVEC array. If SINGVAL is present, the Principal Component is normalized by XSINGVAL on output. This optional argument is useful only if XEIGVEC contains an eigenvector derived from an EOF analysis.

Further Details

The subroutine computes the Principal Component with only one pass through the data by regressing the observations onto the eigenvector of the correlation or covariance matrix of X.

If for some observations in X there is no available data, the corresponding element in XPC is filled with XMISS.

subroutine comp_pc_miss ( x, xeigvec, xpc, xmiss, dimvar, xmean, xstd, xsingval, tol, min_norm )

Purpose

COMP_PC_MISS computes estimates of Principal Components (PC) from a data matrix XX and a set of eigenvectors or singular vectors derived from EOF or MCA analysis when XX contains missing values.

Arguments

X (INPUT) real(stnd), dimension(:,:)
On entry, input submatrix containing size(X,3-DIMVAR) observations on size(X,DIMVAR) variables from the matrix of data XX for which Principal Components are desired. By default, DIMVAR is equal to 1. See description of optional DIMVAR argument for details. If all the data are available at once, X can be the full data matrix XX.
XEIGVEC (INPUT) real(stnd), dimension(:,:)

On entry, XEIGVEC contains selected eigenvectors of the variance-covariance (or correlation) matrix from the data matrix XX or selected left (or right) singular vectors of the SVD of the covariance matrix between the data matrix XX and another data matrix YY.

The shape of XEIGVEC must verify: size( XEIGVEC, 1 ) = size( X, DIMVAR ).

XPC (OUTPUT) real(stnd), dimension(:,:)

On exit, XPC contains the Principal Components derived from X and XEIGVEC.

The shape of XPC must verify:

  • size( XPC, 1 ) = size( X, 3-DIMVAR ) ,
  • size( XPC, 2 ) = size( XEIGVEC, 2 ) .
XMISS (INPUT) real(stnd)
On entry, the missing value indicator. Any value in X which is equal to XMISS is assumed to be missing or invalid. The Principal Components are computed from all the observations where X X are not missing by regression (see Further Details).
DIMVAR (INPUT, OPTIONAL) integer(i4b)

On entry, if DIMVAR is present, DIMVAR is used as follows, if:

  • DIMVAR = 1, the input submatrix X contains size(X,2) observations on size(X,1) variables.
  • DIMVAR = 2, the input submatrix X contains size(X,1) observations on size(X,2) variables.

The default is DIMVAR = 1.

XMEAN (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XMEAN is present, XMEAN contains the variable means and the Principal Components are computed from the centered data matrix X.

The size of XMEAN must verify: size( XMEAN ) = size( X, DIMVAR ).

XSTD (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, if XSTD is present, XSTD contains the variable standard-deviations and the Principal Components are computed from the normalized data matrix X.

The size of XSTD must verify: size( XSTD ) = size( X, DIMVAR ).

XSINGVAL (INPUT, OPTIONAL) real(stnd), dimension(:)

On entry, XSINGVAL must contain the singular values of the covariance (or correlation) matrix from the data matrix XX associated with the eigenvectors in XEIGVEC array. If SINGVAL is present, the Principal Components are normalized by XSINGVAL on output. This optional argument is useful only if XEIGVEC contains eigenvectors derived from an EOF analysis.

The size of XSINGVAL must verify: size( XSINGVAL ) = size( XEIGVEC, 2 ).

TOL (INPUT, OPTIONAL) real(stnd)

On entry, TOL is used to determine the effective rank of the coefficient matrix A for each regression problem, which is then defined as the order of the largest leading triangular submatrix R11 in the QR factorization (with pivoting) of A whose estimated condition number, in the 1-norm, is less than 1/TOL. TOl must be set to the relative precision of the elements in A and B. If each element is correct to, say, 5 digits then TOL=0.00001 should be used.

TOL must not be greater or equal to 1 or less or equal than 0, otherwise the numerical rank of A is determined and the calculations to determine the condition number are not performed. If TOL is absent, the numerical rank of A is determined for each regression problem.

MIN_NORM (INPUT, OPTIONAL) logical(lgl)
On entry, If MIN_NORM=true, minimun 2-norm solutions are computed. If MIN_NORM=false or if MIN_NORM is absent, solutions are computed such that if the j-th column of XEIGVEC is omitted from the basis for the ith observation (regression problem), X[i,j] is set to zero.

Further Details

The subroutine computes the Principal Components with only one pass through the data by regressing the observations onto the eigenvectors of the correlation or covariance matrix of X.

If for some observations in X there is no available data, the corresponding line in XPC is filled with XMISS.

Flag Counter