;----------------------------------------------------------------------------
; $Id: varimax_k58.pro,v 1.3 2002/10/11 21:27:10 johnny Exp johnny $
;+
; NAME:
; VARIMAX_K58
;
; PURPOSE:
; Given an input array of factor loadings (empirical orthogonal fun-
; ctions, i.e. EOFs), rotate the factors using the Varimax criteria
; of Kaiser (1958, 1959).
;
; CATEGORY:
; Statistics.
;
; CALLING SEQUENCE:
; Result = VARIMAX_K58(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.
; Floating array dimensioned (r,p). Unchanged by function.
;
; INPUT KEYWORD PARAMETERS:
; None.
;
; OUTPUT KEYWORD PARAMETERS:
; H: Square root of communality, floating array dimensioned
; (p), of each location (variable). Created/overwritten
; by function.
;
; FRAC_VAR: Fraction of variance accounted for by each r rotated
; EOF (factor), floating array dimensioned (r), assuming
; that the locations (variables) were standardized (see
; p. 203, 213, and other sections in RJ93 for details).
; Created/overwritten by function.
;
; OUTPUTS:
; Result: Column vectors of Result are the first r rotated EOFs
; (factor loadings), floating array dimensioned (r,p).
; Created/overwritten by function.
;
; FILE DEVICE I/O:
; None.
;
; COMMON BLOCKS:
; None.
;
; EXAMPLE:
; This example is taken from Table 7.I in RJ93. The small differences
; in output with Table 7.III are likely due to the use of a different
; criteria in this function for determining when the Varimax solution
; has converged.
;
; Given input factors, conduct a Varimax 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] ]
; B = VARIMAX_K58(A, H=H, Frac_Var=frac_var)
; IDL returns for B:
; -0.0204423 0.938674 -0.340334
; 0.983662 0.0730206 0.134997
; 0.0826106 0.435975 -0.893379
; 0.939901 -0.0965213 -0.309596
; For H^2:
; 0.997354 0.991148 0.995024 0.988580
; For fraction of variance:
; 0.464562 0.271458 0.257007
;
; REFERENCES:
; - Kaiser, H. F. (1958), "The varimax criterion for analytic rotation
; in factor analysis," Psychometrika, Vol. 23, pp. 187-200. Abbrev-
; iated K58 in this function.
; - Kaiser, H. F. (1959), "Computer program for varimax rotation in
; factor analysis," Educ. Psych. Meas., Vol. 19, No. 3, pp. 413-420.
; Abbreviated K59 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.
;
; MODIFICATION HISTORY:
; - 11 Oct 2002: Orig. ver. Johnny Lin, CIRES/University of Colorado.
; Email: air_jlin@yahoo.com. Passed passably adequate tests.
;
; NOTES:
; - Written for IDL 5.3.
; - All keyword parameters are optional unless otherwise stated.
; - No procedures called with _Extra keyword are invoked.
; - Only built-in functions and procedures are called.
;-
; Copyright (c) 2002 Johnny Lin. For licensing, distribution conditions,
; and contact information, see http://www.johnny-lin.com/lib.html.
;----------------------------------------------------------------------------
FUNCTION VARIMAX_K58, in_A $
, H = out_H $
, FRAC_VAR = out_frac_var $
, _EXTRA = extra
; -------------------- Error Check and Parameter Setting --------------------
COMPILE_OPT IDL2
ON_ERROR, 0
epsilon = (MACHAR()).eps ;- set machine precision tolerance
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)
if (dims[0] ne 2) then MESSAGE, 'error--bad input dimensions'
if (ABS(r) gt p) then MESSAGE, 'error--bad r'
; --------------------------- Do Varimax Rotation ---------------------------
;
; Section input: A (undefined upon section output)
; Section output: B, H, H_B
;
; Selected variable key:
; A Unrotated EOFs, dimensioned (r,p). EOFs (factor loadings)
; of input 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 Square root of communality, dimensioned (p)
; H_B Same as H, but as an array in which each element corres-
; ponds to the elements of B and is dimensioned (r,p).
;
; Note: Some of the K59 local variables have the same names as variables
; in this function (e.g. A). To prevent a mix-up, all the K59 inter-
; mediate variables that are one letter long (and not capital X or Y) are
; named in this function as "duplicated" characters, e.g. xx is the same
; as x in K59, AA is the same as A in K59, etc.
;
; Algorithm: Rotation is done 2 EOFs at a time until there is a conver-
; gence of solution (i.e. a set of rotated EOFs B in which the Varimax
; criteria is maximized). Following the idea in K59, convergence of
; solution is assumed when we have tried to rotate through all r*(r-1)/2
; pairs of EOFs in the matrix but all rotation attempts have failed as
; the angle for all pairs is essentially zero. Then the flag not_con-
; verged is set to 0 (false) and the rotation cycle while loop ends.
;--> some initial settings:
B = TEMPORARY(A) ;- set initial B to A
not_converged = 1 ;- init. flag to Varimax solution has not yet converged
;--> calculate square root of communality of each location (variable):
; this is the square root of the sum of the squares of each row of B
; (i.e. sum of squares of loadings for each p location) since the A
; factors (which B currently equals) are uncorrelated (see RJ93, eqns.
; 3.14, 3.1):
H = SQRT( TOTAL(B^2, 1) )
H_B = REBIN(REFORM(H, 1,p), r,p)
;--> calc. B normalized by square root of communality (i.e. b_{ij}/h_{i},
; e.g. in eqn. 7.3 in RJ93) and set as array bh, dimensioned (r,p):
bh = B / H_B
;--> iterate until solution has converged (loop begin):
while not_converged do begin
;- set/reset counter of number of r*(r-1)/2 pairs of EOFs (num_pairs):
num_pairs = r * (r-1) / 2
;- rotate all factors (EOFs) in B two factors at a time: each factor index
; i [0:r-1-1] is paired with each of the remaining factors (index j):
for i=0,r-1-1 do begin ;+ go through factor index 0:r-1-1 (begin)
for j=i+1,r-1 do begin ;+ pair i to "rest" of factors (begin)
;> select pair of factors at factors i,j (notation closely follows K59):
xx = bh[i,*]
yy = bh[j,*]
;> calculate angle of rotation (phi) that maximizes the Varimax function
; for factors i,j (notation closely follows K59; n in K59 is p here;
; see also p. 340 in RJ93); the algorithm used here to calculate phi
; already accounts for maximizing the Varimax function:
uu = xx^2 - yy^2
vv = 2.0 * xx * yy
AA = TOTAL(uu)
BB = TOTAL(vv)
CC = TOTAL(uu^2 - vv^2)
DD = TOTAL(2.0 * uu * vv)
num = DD - (2.0 * AA * BB / p)
den = CC - ( (AA^2 - BB^2) / p )
phi = ATAN(num, den) / 4.0
;> if phi is non-zero, rotate to get new factors (xx_rot and yy_rot,
; which are X and Y in K59), insert into bh; if phi is zero, decrement
; num_pairs by 1 and stop rotation process after one pass through all
; EOF pairs in bh have been shown to have phi=0 rotation angles:
if (ABS(phi) ge 2.0*epsilon) then begin
xx_rot = ( COS(phi) * xx) + (SIN(phi) * yy)
yy_rot = (-SIN(phi) * xx) + (COS(phi) * yy)
bh[i,*] = TEMPORARY(xx_rot)
bh[j,*] = TEMPORARY(yy_rot)
endif else begin
num_pairs = num_pairs - 1
if (num_pairs eq 0) then not_converged = 0
endelse
endfor ;+ pair i to "rest" of factors (end)
endfor ;+ go through factor index 0:r-1-1 (end)
;--> iterate until solution has converged (loop end):
endwhile
;--> denormalize (and deallocate) bh to get final B by multiplying bh by
; the square root of communality:
B = TEMPORARY(bh) * H_B
; ---------------------- Calculate Fraction of Variance ---------------------
;
; Section input: B
; Section output: fvar
;
; Algorithm: Calculate fraction of variance (fvar) for each EOF/PC
; (factor) for the B matrix, as described on the bottom of p. 213,
; Table 7.I-7.III, and top of p. 215 of RJ93.
fvar = TOTAL(B^2, 2) / FLOAT(p)
; ------------------------------ Prepare Output -----------------------------
out_H = TEMPORARY(H) ;- output sqrt of communality (optional)
out_frac_var = TEMPORARY(fvar) ;- output frac. of var. of EOF/PCs (optional)
Result = TEMPORARY(B) ;- output rotated EOFs
RETURN, Result
END ;=== end of function ===
; ========== end file ==========