# SciRuby/statsample

lib/statsample/factor/rotation.rb

### Summary

B
6 hrs
###### Test Coverage
``````module Statsample
module Factor
# Base class for component matrix rotation.
#
# == Reference:
# * SPSS Manual
# * Lin, J. (2007). VARIMAX_K58 [Source code]. [http://www.johnny-lin.com/idl_code/varimax_k58.pro]
#
# Use subclasses Varimax, Equimax or Quartimax for desired type of rotation
#   Use:
#   a = Matrix[ [ 0.4320,  0.8129,  0.3872]
#     , [ 0.7950, -0.5416,  0.2565]
#     , [ 0.5944,  0.7234, -0.3441]
#     , [ 0.8945, -0.3921, -0.1863] ]
#   rotation = Statsample::Factor::Varimax(a)
#   rotation.iterate
#   p rotation.rotated
#   p rotation.component_transformation_matrix
#
class Rotation
EPSILON=1e-15
MAX_ITERATIONS=25
include Summarizable
include DirtyMemoize
# Maximum number of iterations
attr_accessor :max_iterations
# Maximum precision
attr_accessor :epsilon
attr_accessor :use_gsl
dirty_writer :max_iterations, :epsilon
dirty_memoize :iterations, :rotated, :component_transformation_matrix, :h2

def initialize(matrix, opts=Hash.new)
@name=_("%s rotation") % rotation_name
@matrix=matrix
@n=@matrix.row_size # Variables, p on original
@m=@matrix.column_size # Factors, r on original
@component_transformation_matrix=nil
@max_iterations=MAX_ITERATIONS
@epsilon=EPSILON
@rotated=nil
@h2=(@matrix.collect {|c| c**2} * Matrix.column_vector([1]*@m)).column(0).to_a
@use_gsl=Statsample.has_gsl?
opts.each{|k,v|
self.send("#{k}=",v) if self.respond_to? k
}
end
def report_building(g)
g.section(:name=>@name) do |s|
s.parse_element(rotated)
s.parse_element(component_transformation_matrix)
end
end
alias_method :communalities, :h2
alias_method :rotated_component_matrix, :rotated
def compute
iterate
end
# Start iteration
def iterate
k_matrix=@use_gsl ? GSL::Matrix : ::Matrix
t=k_matrix.identity(@m)
b=(@use_gsl ? @matrix.to_gsl : @matrix.dup)
h=k_matrix.diagonal(*@h2).collect {|c| Math::sqrt(c)}
h_inverse=h.collect {|c| c!=0 ? 1/c : 0 }
bh=h_inverse * b
@not_converged=true
@iterations=0
while @not_converged
break if @iterations>@max_iterations
@iterations+=1
#puts "Iteration #{iterations}"
num_pairs=@m*(@m-1).quo(2)
(0..(@m-2)).each do |i| #+ go through factor index 0:r-1-1 (begin)
((i+1)..(@m-1)).each do |j| #+ pair i to "rest" of factors (begin)

xx = bh.column(i)
yy = bh.column(j)
tx = t.column(i)
ty = t.column(j)

uu = @n.times.collect {|var_i| xx[var_i]**2-yy[var_i]**2}
vv = @n.times.collect {|var_i| 2*xx[var_i]*yy[var_i]}

a  = @n.times.inject(0) {|ac,var_i| ac+ uu[var_i] }
b  = @n.times.inject(0) {|ac,var_i| ac+ vv[var_i] }
c  = @n.times.inject(0) {|ac,var_i| ac+ (uu[var_i]**2 - vv[var_i]**2) }
d  = @n.times.inject(0) {|ac,var_i| ac+ (2*uu[var_i]*vv[var_i]) }
num=x(a,b,c,d)
den=y(a,b,c,d)
phi=Math::atan2(num,den) / 4.0
# puts "#{i}-#{j}: #{phi}"

if(Math::sin(phi.abs) >= @epsilon)
xx_rot=( Math::cos(phi)*xx)+(Math::sin(phi)*yy)
yy_rot=((-Math::sin(phi))*xx)+(Math::cos(phi)*yy)

tx_rot=( Math::cos(phi)*tx)+(Math::sin(phi)*ty)
ty_rot=((-Math::sin(phi))*tx)+(Math::cos(phi)*ty)

bh=bh.to_a

@n.times {|row_i|
bh[row_i][i] = xx_rot[row_i]
bh[row_i][j] = yy_rot[row_i]
}
t=t.to_a
@m.times {|row_i|
t[row_i][i]=tx_rot[row_i]
t[row_i][j]=ty_rot[row_i]
}
#if @use_gsl
bh=k_matrix.[](*bh)
t=k_matrix.[](*t)
#else
#  bh=Matrix.rows(bh)
#  t=Matrix.rows(t)

#end
else
num_pairs=num_pairs-1
@not_converged=false if num_pairs==0
end # if
end #j
end #i
end # while
@rotated=h*bh
@rotated.extend CovariateMatrix
@rotated.name=_("Rotated Component matrix")

if @matrix.respond_to? :fields_x
@rotated.fields_x = @matrix.fields_x
else
@rotated.fields_x = @n.times.map {|i| "var_#{i+1}"}
end
if @matrix.respond_to? :fields_y
@rotated.fields_y = @matrix.fields_y
else
@rotated.fields_y = @m.times.map {|i| "var_#{i+1}"}
end

@component_transformation_matrix=t
@component_transformation_matrix.extend CovariateMatrix
@component_transformation_matrix.name=_("Component transformation matrix")

if @matrix.respond_to? :fields_y
@component_transformation_matrix.fields = @matrix.fields_y

else
@component_transformation_matrix.fields = @m.times.map {|i| "var_#{i+1}"}
end

@rotated
end

end
class Varimax < Rotation
def x(a,b,c,d)
d-(2*a*b / @n.to_f)
end
def y(a,b,c,d)
c-((a**2-b**2) / @n.to_f)
end
def rotation_name
"Varimax"
end
end
class Equimax < Rotation
def x(a,b,c,d)
d-(@m*a*b / @n.to_f)
end
def y(a,b,c,d)
c-@m*((a**2-b**2) / (2*@n.to_f))
end
def rotation_name
"Equimax"
end

end
class Quartimax < Rotation
def x(a,b,c,d)
d
end
def y(a,b,c,d)
c
end
def rotation_name
"Quartimax"
end

end
end
end``````