1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
|
#Varimax matrix rotation
#Coded in Python by Sam Webb in 2003
#Adapted from MATLAB code modified by Tevor Park in 2002 from original code by J.O. Ramsay
#
from Numeric import *
from math import *
import MLab
import LinearAlgebra
def angle(z):
a=z.real
b=z.imag
return arctan2(b,a)
def var(z):
s=MLab.std(z)
return s**2
def vcomplex(a,b):
return a+b*1j
def varimax(amat,target_basis=0):
MAX_ITER=50
EPSILON=1e-10
amatd=amat.shape
#make sure the problem is 2D
if len(amatd)!=2:
raise ValueError, 'AMAT must be 2-dimensional'
n=amatd[0]
k=amatd[1]
rotm=identity(k)*1.0
#if singular
if k==1:
return
#error check target_basis
if target_basis!=0:
if len(target_basis.shape)!=2:
raise ValueError, 'TARGET_BASIS must be 2-dimensional'
if alltrue(target_basis.shape==(n,n)):
amat=LinearAlgebra.solve_linear_equations(target_basis,amat)
else:
raise ValueError, 'TARGET_BASIS must be a basis for the column space'
else:
target_basis=identity(n)*1.0
#on to the guts
varnow=sum(var(amat**2))
not_converged=1
iter=0
while not_converged and iter<MAX_ITER:
for j in range(0,k-1):
for l in range(j+1,k):
#find optimal planar rotation angle for column j,l
#break expression into parts
c1=vcomplex(amat[:,j],amat[:,l])
c1=n*sum(c1**4)
c2=vcomplex(amat[:,j],amat[:,l])
c2=sum(c2**2)
phi_max=angle(c1-c2**2)/4
sub_rot=array([[cos(phi_max),-sin(phi_max)],[sin(phi_max),cos(phi_max)]])
atemp=take(amat,(j,l),1)
rtemp=take(rotm,(j,l),1)
atemp=matrixmultiply(atemp,sub_rot)
rtemp=matrixmultiply(rtemp,sub_rot)
put(amat,range(j,n*k,k),atemp[:,0])
put(amat,range(l,n*k,k),atemp[:,1])
put(rotm,range(j,k*k,k),rtemp[:,0])
put(rotm,range(l,k*k,k),rtemp[:,1])
varold=varnow
varnow=sum(var(amat**2))
if varnow==0:
return
not_converged=((varnow-varold)/varnow > EPSILON)
iter=iter+1
if iter>=MAX_ITER:
print 'WARNING: Maximum number of iterations exceeded in varimax rotation'
opt_amat=matrixmultiply(target_basis,amat)
print "Total varimax iterations: "+str(iter)
return rotm,opt_amat
|
