Viewing contents of file '../idllib/idl_5.2/lib/correlate.pro'
;$Id: correlate.pro,v 1.13.4.1 1999/01/16 16:38:07 scottm Exp $
;
; Copyright (c) 1994-1999, Research Systems, Inc. All rights reserved.
; Unauthorized reproduction prohibited.
;+
; NAME:
; CORRELATE
;
; PURPOSE:
; This function computes the linear Pearson correlation coefficient
; of two vectors or the Correlation Matrix of an M x N array.
; Alternatively, this function computes the covariance of two vectors
; or the Covariance Matrix of an M x N array.
;
; CATEGORY:
; Statistics.
;
; CALLING SEQUENCE:
; Result = Correlate(X [,Y])
;
; INPUTS:
; X: A vector or an M x N array of type integer, float or double.
;
; Y: A vector of type integer, float or double. If X is an M x N
; array, this parameter should not be supplied.
;
; KEYWORD PARAMETERS:
; COVARIANCE: If set to a non-zero value, the sample covariance is
; computed.
;
; DOUBLE: If set to a non-zero value, computations are done in
; double precision arithmetic.
;
; RESTRICTIONS:
; If X is an M x N array, then Y should not be supplied;
; Result = Correlate(X)
;
; EXAMPLES:
; Define the data vectors.
; x = [65, 63, 67, 64, 68, 62, 70, 66, 68, 67, 69, 71]
; y = [68, 66, 68, 65, 69, 66, 68, 65, 71, 67, 68, 70]
;
; Compute the linear Pearson correlation coefficient of x and y.
; result = correlate(x, y)
; The result should be 0.702652
;
; Compute the covariance of x and y.
; result = correlate(x, y, /covariance)
; The result should be 3.66667
;
; Define an array with x and y as its columns.
; a = [transpose(x), transpose(y)]
; Compute the correlation matrix.
; result = correlate(a)
; The result should be [[1.000000, 0.702652]
; [0.702652, 1.000000]]
;
; Compute the covariance matrix.
; result = correlate(a, /covariance)
; The result should be [[7.69697, 3.66667]
; [3.66667, 3.53788]]
;
; PROCEDURE:
; CORRELATE computes the linear Pearson correlation coefficient of
; two vectors. If the vectors are of unequal length, the longer vector
; is truncated to the length of the shorter vector. If X is an M x N
; array, M-columns by N-rows, the result will be an M x M array of
; linear Pearson correlation coefficients with the iTH column and jTH
; row element corresponding to the correlation of the iTH and jTH
; columns of the M x N array. The M x M array of linear Pearson
; correlation coefficients (also known as the Correlation Matrix) is
; always symmetric and contains 1s on the main diagonal. The Covariance
; Matrix is also symmetric, but is not restricted to 1s on the main
; diagonal.
;
; REFERENCE:
; APPLIED STATISTICS (third edition)
; J. Neter, W. Wasserman, G.A. Whitmore
; ISBN 0-205-10328-6
;
; MODIFICATION HISTORY:
; Written by: DMS, RSI, Sept 1983
; Modified: GGS, RSI, July 1994
; Added COVARIANCE keyword.
; Included support for matrix input.
; New documentation header.
; Modified: GGS, RSI, April 1996
; Included support for scalar and unequal length vector
; inputs. Added checking for undefined correlations and
; support of IEEE NaN (Not a Number).
; Added DOUBLE keyword.
; Modified keyword checking and use of double precision.
;-
FUNCTION idl_corr_total, Arg, Double = Double
Type = SIZE(Arg)
if N_ELEMENTS(Double) eq 0 then $
Double = (Type[Type[0]+1] eq 5)
ArgSum = TOTAL(Arg, Double = Double)
if Type[Type[0]+1] eq 5 and Double eq 0 then $
RETURN, FLOAT(ArgSum) else RETURN, ArgSum
END
FUNCTION Cov_Mtrx, X, Double = Double
Nx = SIZE(x)
if N_ELEMENTS(Double) eq 0 then Double = (Nx[Nx[0]+1] eq 5)
if Double eq 0 then one = 1.0 else one = 1.0d
if nx[0] le 1 then RETURN, one
VarXi = X - TOTAL(x,2)/nx[2] # REPLICATE(one, nx[2])
VarXi = (VarXi # TRANSPOSE(VarXi)) / (nx[2]-1)
if Double eq 0 then RETURN, FLOAT(VarXi) else RETURN, VarXi
end
FUNCTION CRR_MTRX, X, Double = Double
Nx = SIZE(x)
if N_ELEMENTS(Double) eq 0 then Double = (Nx[Nx[0]+1] eq 5)
if Double eq 0 then one = 1.0 else one = 1.0d
if Nx[0] eq 0 or Nx[0] eq 1 then RETURN, one
VarXi = x - TOTAL(x,2)/nx[2] # REPLICATE(one, Nx[2])
SS = VarXi^2 # REPLICATE(one, Nx[2])
SS = SS # SS
Tol = sqrt((machar(DOUBLE=double)).xmin)
iSS = WHERE(ABS(SS) le Tol, nSS) ;Zero denominator signals undefined
;Correlation.
if nSS ne 0 then begin
SS[iSS] = one
cm = (VarXi # TRANSPOSE(VarXi))/SQRT(SS)
cm[iSS] = !VALUES.F_NAN
if Double eq 0 then RETURN, FLOAT(cm) else RETURN, cm
endif else $
if Double eq 0 then RETURN, FLOAT((VarXi # TRANSPOSE(VarXi))/SQRT(SS)) $
else RETURN, (VarXi # TRANSPOSE(VarXi))/SQRT(SS)
end
FUNCTION Correlate, X, Y, Covariance = Covariance, Double = Double
ON_ERROR, 2 ;Return to caller if an error occurs.
if N_PARAMS() eq 2 then begin ;Vector inputs.
Sx = SIZE(x) & Sy = SIZE(y)
if N_ELEMENTS(Double) eq 0 then $
Double = (Sx[Sx[0]+1] eq 5) or (Sy[Sy[0]+1] eq 5)
Nx = Sx[Sx[0]+2]
Ny = Sy[Sy[0]+2]
;If not of equal length, take shortest.
if Nx lt Ny then begin
yMem = Y[Nx:*] & Y = Y[0:Nx-1]
endif else $
if Ny lt Nx then begin
xMem = X[Ny:*] & X = X[0:Ny-1]
endif
;Means.
sLength = MIN([Nx, Ny])
xMean = idl_corr_total(X, Double = Double) / sLength
yMean = idl_corr_total(Y, Double = Double) / sLength
;Deviations.
xDev = X - xMean
yDev = Y - yMean
;If not of equal length, restore input vector.
if Nx lt Ny then Y = [Y, yMem] else $
if Ny lt Nx then X = [X, xMem]
if KEYWORD_SET(Covariance) eq 0 then begin ;Correlation.
Denom = idl_corr_total(xDev^2, Double = Double) * $
idl_corr_total(yDev^2, Double = Double)
Tol = sqrt((machar(DOUBLE=double)).xmin)
iDenom = WHERE(ABS(Denom) le Tol, nDenom) ;Zero denominator signals
;undefined Correlation.
if nDenom ne 0 then begin
Denom[iDenom] = 1.0
cr = idl_corr_total(xDev * yDev, Double = Double)/SQRT(Denom)
cr[iDenom] = !VALUES.F_NAN
RETURN, cr
endif else RETURN, idl_corr_total(xDev * yDev, Double = Double) / $
SQRT(Denom)
endif else begin ;Covariance.
if sLength eq 1 then RETURN, !VALUES.F_NAN $
else RETURN, idl_corr_total(xDev * yDev, Double = Double) / (sLength-1)
endelse
endif else begin ;Array input.
if KEYWORD_SET(Covariance) eq 0 then RETURN, $
CRR_MTRX(X, Double = Double) $ ;Correlation Matrix.
else RETURN, COV_MTRX(X, Double = Double) ;Covariance Matrix.
endelse
END