;----------------------------------------------------------------------------
; $Id: promax_hw64.pro,v 1.10 2002/11/06 02:59:39 johnny Exp $
;+
; NAME:
; PROMAX_HW64
;
; PURPOSE:
; Given an input array of factor loadings (empirical orthogonal fun-
; ctions, i.e. EOFs), rotate the factors using a method almost iden-
; tical to the Promax method of Hendrickson and White (1964).
;
; CATEGORY:
; Statistics.
;
; CALLING SEQUENCE:
; Result = PROMAX_HW64(in_A)
;
; INPUTS:
; in_A: Column vectors of in_A are the first r unrotated EOFs
; (factor loadings) of the input data (of n observations
; at p locations) from which these EOFs are derived.
; One can generate in_a by taking the transpose of the
; values returned in the Coefficients keyword in the built-
; in function PCOMP. Floating array dimensioned (r,p).
; Unchanged by function.
;
; INPUT KEYWORD PARAMETERS:
; M: Power m to raise "target matrix" to (while retaining sign
; of said target matrix) to create the "ideal pattern
; matrix" (RJ93) which decreases high loadings less than
; low loadings. See eqn. 28 in R86. If keyword is not
; defined, m defaults to 2 (following R86, p. 320). Un-
; changed by function.
;
; OUTPUT KEYWORD PARAMETERS:
; P_P: Column vectors of P_P are the first r Promax oblique
; primary factor pattern EOFs (factor loadings), floating
; array dimensioned (r,p). This is P_p in RJ93, p. 221.
; Created/overwritten by function.
;
; P_R: Column vectors of P_R are the first r Promax oblique
; reference pattern EOFs (factor loadings), floating array
; dimensioned (r,p). This is P_r in RJ93, p. 221. Created/
; overwritten by function.
;
; PHI_P: Correlation between primary factors, floating array di-
; mentioned (r,r). Each [i,j] element is the correlation
; between factors i and j. This is phi_p in RJ93, p. 221.
; Created/overwritten by function.
;
; PHI_R: Correlation between reference axes, floating array dimen-
; tioned (r,r). Each [i,j] element is the correlation
; between reference axes i and j. This is phi_r in RJ93,
; p. 221. Created/overwritten by function.
;
; S_P: Column vectors of S_P are the first r Promax oblique
; primary factor structure EOFs (factor loadings), floating
; array dimensioned (r,p). This is S_p in RJ93, p. 221.
; Created/overwritten by function.
;
; S_R: Column vectors of S_R are the first r Promax oblique
; reference structure EOFs (factor loadings), floating
; array dimensioned (r,p). This is S_r in RJ93, p. 220
; and is the same as the matrix passed out via Result.
; Created/overwritten by function.
;
; T_P: Primary structure transformation matrix, floating
; array dimensioned (r,r). This is T_p in RJ93, p. 221.
; Created/overwritten by function.
;
; T_R: Reference structure transformation matrix, floating
; array dimensioned (r,r). This is T_r in RJ93, p. 220.
; Created/overwritten by function.
;
; OUTPUTS:
; Result: Column vectors of Result are the first r Promax oblique
; reference structure EOFs (factor loadings), floating
; array dimensioned (r,p). This is S_r in RJ93, p. 220
; and is the same as the matrix passed out by output
; keyword S_R. Created/overwritten by function.
;
; FILE DEVICE I/O:
; None.
;
; COMMON BLOCKS:
; None.
;
; EXAMPLE:
; This example uses input data from Table 7.I in RJ93 and gives
; results that are similar to Tables 7.IV and 7.V. The differences
; in output with those tables are probably because the Varimax routine
; used in this function probably has different iteration criteria than
; RJ93, and it appears in calculating Tables 7.IV and 7.V the B matrix
; was not forced row normalized in RJ93 (pp. 220, 343, and see note
; below in this header block).
;
; Given input factors, conduct a Promax rotation:
; A = [ [ 0.4320, 0.8129, 0.3872] $
; , [ 0.7950, -0.5416, 0.2565] $
; , [ 0.5944, 0.7234, -0.3441] $
; , [ 0.8945, -0.3921, -0.1863] ]
; S_r = PROMAX_HW64(A, M=4, P_P=P_p, S_P=S_p)
; IDL returns for P_p:
; 0.0108750 0.912841 -0.160862
; 1.02958 0.162965 0.254964
; -0.0149609 0.168635 -0.910244
; 0.902053 -0.182991 -0.301612
; For S_p:
; -0.00528595 0.989576 -0.601337
; 0.975345 -0.0105829 -0.00968674
; 0.141439 0.606637 -0.988544
; 0.965725 -0.0828474 -0.376990
;
; REFERENCES:
; - Hendrickson, A. E., and P. O. White (1964), "Promax: A quick
; method for rotation to oblique simple structure," British J.
; Statist. Psychol., Vol. 17, pp. 65-70. Abbreviated HW64 in this
; function.
; - Reyment, Richard A., and K. G. Joreskog (1993), Applied Factor
; Analysis in the Natural Sciences, 2nd Edition (Cambridge: Cambridge
; University Press), 371 pp. Abbreviated RJ93 in this function. See
; also MATLAB code by L. F. Marcus in the Appendix of this book.
; - Richman, M. B. (1986), "Rotation of principal components," Journal
; of Climatology, Vol. 6, pp. 293-335. Abbreviated R86 in this
; function.
;
; MODIFICATION HISTORY:
; - 18 Oct 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado.
; Email: air_jlin@yahoo.com. Passed passably adequate tests.
; - 4 Nov 2002: Fix variance conservation problem. Passed visual
; and reasonable values tests.
;
; NOTES:
; - Written for IDL 5.3.
; - The output structure and pattern matrices are all based upon a
; row normalized target matrix B and column normalized reference
; and pattern transformation matrices. Thus, even if the input
; EOFs are non-normalized, the output EOFs will be "normalized"
; (see references for this function).
; - The form of the HW64's Promax method used is that as described in
; RJ93 and R86. The "target matrix" is chosen to be the Varimax ro-
; tated EOFs (factor loadings) matrix.
; - In R86, p. 320, the columns and rows of B (called A* in R86) are
; normalized. In testing RPCA_LA00, which uses PROMAX_HW64, it was
; found that normalizing columns of B causes total variance to not be
; conserved during rotation (whether one uses the method of R86 or
; regular normalization). Thus, in this implementation of Promax,
; only rows are normalized in B; this is the major deviation from
; the canonical Promax method.
; - RJ93, pp. 215-220 has a detailed discussion of what "primary" and
; "reference" factors are as well as "structure" and "pattern" loadings.
; - This function was tested using RJ93, Table 7.III as values for the
; Varimax target matrix (instead of calculating the Varimax matrix using
; VARIMAX_K58), and without row-column normalizing B at the beginning
; of the Promax procedure. This yielded results nearly identical to
; Tables 7.IV and 7.V, which is the basis of the explanation of the
; example discrepancy with RJ93 discussed above.
; - All keyword parameters are optional unless otherwise stated.
; - No procedures called with _Extra keyword are invoked.
; - Non-built-in functions and procedures called: VARIMAX_K58.
;-
; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions,
; and contact information, see http://www.johnny-lin.com/lib.html.
;----------------------------------------------------------------------------
FUNCTION PROMAX_HW64, in_A $
, M = in_m $
, P_P = out_P_p $
, P_R = out_P_r $
, PHI_P = out_phi_p $
, PHI_R = out_phi_r $
, S_P = out_S_p $
, S_R = out_S_r $
, T_P = out_T_p $
, T_R = out_T_r $
, _EXTRA = extra
; -------------------- Error Check and Parameter Setting --------------------
COMPILE_OPT IDL2
ON_ERROR, 0
A = in_A ;- protect input data
dims = SIZE(A) ;- dimensions of input data
r = dims[1] ;- set number of retained EOFs (factor loadings)
; (in RJ93 eqn. 7.2 etc., this is called k):
p = dims[2] ;- set number of locations (in RJ93 these are called
; variables)
in_type = dims[3] ;- set input data type
if (dims[0] ne 2) then MESSAGE, 'error--bad input dimensions'
if (ABS(r) gt p) then MESSAGE, 'error--bad r'
if (N_ELEMENTS(in_m) eq 0) then begin ;- set value of exponent m (see
m = 2 ; descrip. of keyword M; default
endif else begin ; for m is 2
m = in_m
endelse
; ---------- Do Varimax Rotation To Get Preliminary Promax Solution ---------
;
; Section input: A (undefined on section exit)
; Section output: B, H_B
;
; Selected variable key:
; A First r unrotated EOFs (factor loadings), dimensioned (r,p).
; EOFs (factor loadings) of orig. data are column vectors of A.
; B Rotated EOFs, dimensioned (r,p). The first r rotated EOFs
; (factor loadings) are column vectors of B.
; H_B The square root of communality as an array in which each
; element corresponds to the elements of B. Dimensioned
; (r,p).
B = VARIMAX_K58(TEMPORARY(A), H=H)
H_B = REBIN(REFORM(H, 1,p), r,p)
; -------------------------- Find Promax Solution ---------------------------
;
; Section input: B, H_B
; Section output: B (changed), P_p, P_r, S_p, S_r, T_p, T_r
;
; The obliquely rotated EOFs are column vectors of the output (P_p, P_r,
; etc.)
;
; Selected variable key:
; B Target matrix (RJ93, p. 220) which are the Varimax rotated
; EOFs, dimensioned (r,p). On section input, the first r
; rotated EOFs (factor loadings) are column vectors of B and
; called A* in R86, p. 320. Rows forced to be normalized in
; this section.
; B_star "Ideal pattern matrix" (RJ93, p. 220), dimensioned (r,p).
; Calculating using method described in R86, p. 320 (R86
; calls B_star "B"). Note that RJ93 (p. 220) incorrectly
; implies B_star=B^m, without retaining the sign of B (R86
; is correct, if one compares back to HW64).
; inv_stat Flag storing whether any of the INVERT functions in the
; section did not work correctly. If non-zero then fctn.
; will stop executing.
; max_col Maximum of abs. abs. value in each col. of B (row normal-
; ized). Dimensioned (r,p) where each element corresponds
; one-to-one with B.
; P_p See description of keyword P_P.
; P_r See description of keyword P_R.
; S_p See description of keyword S_P.
; S_r See description of keyword S_R.
; T_p See description of keyword T_P.
; T_p_trans Transpose of T_p.
; T_r See description of keyword T_R.
; T_r_trans Transpose of T_r.
;
; Algorithm: Most notation is similar to RJ93, p. 220. Algorithms are
; from both RJ93 and R86; see the sub-sections for details. See the
; "Notes" section in the top header for additional information re the
; Promax terminology.
;- set initial inv_stat to 0:
inv_stat = 0
;- normalization of rows of B by H_B (B is called A* in R86):
B = TEMPORARY(B) / H_B
;- calculate "ideal pattern matrix" (B_star in RJ93, p. 220; algorithm
; follows R86 eqn. 28):
neg_pts = WHERE(B lt 0.0, cnt_neg) ;+ indices of which elements are
pos_pts = WHERE(B ge 0.0, cnt_pos) ; what sign in the B matrix
B_star = B^m
if (cnt_neg gt 0) then B_star[neg_pts] = (ABS(B_star))[neg_pts] * (-1.0)
if (cnt_pos gt 0) then B_star[pos_pts] = (ABS(B_star))[pos_pts]
;- calculate transformation matrix T_r (RJ93, p. 220) and normalize col.
; (RJ93, p. 343); also calculate T_r_trans:
T_r = INVERT(TRANSPOSE(B) ## B, status) ## TRANSPOSE(B) ## B_star
inv_stat = inv_stat > status
tmp = MAKE_ARRAY(Dimension=SIZE(T_r, /Dimensions), Type=in_type)
idx = LINDGEN((SIZE(T_r, /Dimensions))[0])
tmp[idx,idx] = (TRANSPOSE(T_r) ## T_r)[idx,idx]
T_r = TEMPORARY(T_r) ## SQRT( INVERT(TEMPORARY(tmp), status) )
inv_stat = inv_stat > status
T_r_trans = TRANSPOSE(T_r)
;- calculate transformation matrix T_p (RJ93, p. 220) with the rows of
; transpose of T_p (T_p_trans) normalized (RJ93, p. 343):
T_p_trans = INVERT(T_r, status)
inv_stat = inv_stat > status
tmp = MAKE_ARRAY(Dimension=SIZE(T_p_trans, /Dimensions), Type=in_type)
idx = LINDGEN((SIZE(T_p_trans, /Dimensions))[0])
tmp[idx,idx] = (T_p_trans ## TRANSPOSE(T_p_trans))[idx,idx]
T_p_trans = INVERT(SQRT(TEMPORARY(tmp)), status) ## TEMPORARY(T_p_trans)
inv_stat = inv_stat > status
T_p = TRANSPOSE(T_p_trans)
;- calculate structures (RJ93, p. 221): primary factor structure S_p
; and oblique reference structure S_r, both dimensioned (r,p):
S_p = B ## T_p
S_r = B ## T_r
;- calculate patterns (RJ93, p. 221): primary factor pattern P_p and
; oblique reference pattern P_r, both dimensioned (r,p):
P_p = B ## INVERT(T_p_trans, status)
inv_stat = inv_stat > status
P_r = B ## INVERT(T_r_trans, status)
inv_stat = inv_stat > status
;- check if any of the matrix inverses are bad:
if (inv_stat ne 0) then MESSAGE, 'error--bad inversions'
; ---------- Calc. Correl. Betw. Primary Factors and Reference Axes ---------
;
; Section input: T_p, T_p_trans, T_r, T_r_trans
; Section output: phi_p, phi_r
;
; Selected variable key:
; phi_p Correlations between primary factors, dimensioned (r,r).
; phi_r Correlations between reference axes, dimensioned (r,r).
;
; Algorithm: See RJ93, p. 221 for details.
phi_p = T_p_trans ## T_p
phi_r = T_r_trans ## T_r
; ------------------------------ Prepare Output -----------------------------
out_P_p = TEMPORARY(P_p) ;- keyword P_P (optional)
out_P_r = TEMPORARY(P_r) ;- keyword P_R (optional)
out_phi_p = TEMPORARY(phi_p) ;- keyword PHI_P (optional)
out_phi_r = TEMPORARY(phi_r) ;- keyword PHI_R (optional)
out_S_p = TEMPORARY(S_p) ;- keyword S_P (optional)
out_S_r = TEMPORARY(S_r) ;- keyword S_R (optional)
out_T_p = TEMPORARY(T_p) ;- keyword T_P (optional)
out_T_r = TEMPORARY(T_r) ;- keyword T_R (optional)
Result = out_S_r ;- main function output
RETURN, Result
END ;=== end of function ===
; ========== end file ==========