;$Id: p_correlate.pro,v 1.12 2001/01/15 22:28:08 scottm Exp $
;
; Copyright (c) 1994-2001, Research Systems, Inc.  All rights reserved.
;       Unauthorized reproduction prohibited.
;+
; NAME:
;       P_CORRELATE
;
; PURPOSE:
;       This function computes the partial correlation coefficient of a
;       dependent variable and one particular independent variable when
;       the effects of all other variables involved are removed.
;
; CATEGORY:
;       Statistics.
;
; CALLING SEQUENCE:
;       Result = P_correlate(X, Y, C)
;
; INPUTS:
;       X:    An n-element vector of type integer, float or double that
;             specifies the independent variable data.
;
;       Y:    An n-element vector of type integer, float or double that
;             specifies the dependent variable data.
;
;       C:    An array of type integer, float or double that specifies the
;             independent variable data whose effects are to be removed.
;             The columns of this two dimensional array correspond to the
;             n-element vectors of independent variable data.
;
; KEYWORD PARAMETERS:
;    DOUBLE:  If set to a non-zero value, computations are done in
;             double precision arithmetic.
;
; EXAMPLES:
;       Define the data vectors.
;         x0 = [64, 71, 53, 67, 55, 58, 77, 57, 56, 51, 76, 68]
;         x1 = [57, 59, 49, 62, 51, 50, 55, 48, 52, 42, 61, 57]
;         x2 = [ 8, 10,  6, 11,  8,  7, 10,  9, 10,  6, 12,  9]
;
;       Compute the partial correlation of x0 and x1 with the effects of
;       x2 removed. The result should be 0.533469
;         result = p_correlate(x0, x1, reform(x2, 1, n_elements(x2)))
;
;       Compute the partial correlation of x0 and x2 with the effects of
;       x1 removed. The result should be 0.334572
;         result = p_correlate(x0, x2, reform(x1, 1, n_elements(x1)))
;
;       Compute the partial correlation of x1 and x2 with the effects of
;       x0 removed. The result should be 0.457907
;         result = p_correlate(x1, x2, reform(x0, 1, n_elements(x0)))
;
; REFERENCE:
;       APPLIED STATISTICS (third edition)
;       J. Neter, W. Wasserman, G.A. Whitmore
;       ISBN 0-205-10328-6
;
; MODIFICATION HISTORY:
;       Modified by:  GGS, RSI, July 1994
;                     Minor changes to code. New documentation header.
;       Modified by:  GGS, RSI, August 1996
;                     Added DOUBLE keyword.
;                     Modified keyword checking and use of double precision.
;-

FUNCTION P_Correlate, X, Y, C, Double = Double

  COMPILE_OPT idl2

  ON_ERROR, 2  ;Return to caller if an error occurs.

  Sx = SIZE(x)  &  Sy = SIZE(y)  &  Sc = SIZE(c)

  if Sx[Sx[0]+2] ne Sy[Sy[0]+2] then MESSAGE, $
    "X and Y must have the same number of elements."

  if Sc[0] ne 2 then MESSAGE, $
    "C parameter must be a two-dimensional array."

  ;Check row dimension of C.
  if Sx[Sx[0]+2] ne Sc[Sc[0]] then MESSAGE, $
    "Incompatible arrays."

  if N_ELEMENTS(Double) eq 0 then $
    Double = (Sx[Sx[0]+1] eq 5) or (Sy[Sy[0]+1] eq 5) or (Sc[Sc[0]+1] eq 5)

  if Sc[1] eq 1 then begin
    p = [CORRELATE(X, Y, Double = Double), $
         CORRELATE(X, C, Double = Double), $
         CORRELATE(Y, C, Double = Double)]
    if (p[1] ne 1 and p[2] ne 1) then $
      RETURN, (p[0] - p[1] * p[2])/SQRT((1 - p[1]^2) * (1 - p[2]^2)) $
      else RETURN, 0 * p
  endif else begin
    ;Vector of weights.
    if Double eq 0 then Wts = REPLICATE(1.0, Sc[2]) $
    else Wts = REPLICATE(1.0d, Sc[2])
	dummy = REGRESS(C, Y, DOUBLE=double, MCOR=p0)
	dummy = REGRESS([C, TRANSPOSE(X)], Y, DOUBLE=double, MCOR=p1)
    if Double eq 0 then begin
      p0 = FLOAT(p0)
      p1 = FLOAT(p1)
    endif
    p0 = 1 - p0^2
    p1 = 1 - p1^2
    if p0 eq 0 then $
      RETURN, 0 * p $
      else RETURN, SQRT((p0 - p1)/p0)
  endelse

END