Class: Statsample::Factor::PCA

Inherits:
Object
  • Object
show all
Includes:
Summarizable
Defined in:
lib/statsample/factor/pca.rb

Overview

Principal Component Analysis (PCA) of a covariance or correlation matrix..

NOTE: Sign of second and later eigenvalues could be different using Ruby or GSL, so values for PCs and component matrix should differ, because extendmatrix and gsl’s methods to calculate eigenvectors are different. Using R is worse, cause first eigenvector could have negative values! For Principal Axis Analysis, use Statsample::Factor::PrincipalAxis

Usage:

require 'statsample'
a = Daru::Vector.new([2.5, 0.5, 2.2, 1.9, 3.1, 2.3, 2.0, 1.0, 1.5, 1.1])
b = Daru::Vector.new([2.4,0.7,2.9,2.2,3.0,2.7,1.6,1.1,1.6,0.9])
ds = Daru::DataFrame.new({:a => a,:b => b})
cor_matrix = Statsample::Bivariate.correlation_matrix(ds)
pca=  Statsample::Factor::PCA.new(cor_matrix)
pca.m
=> 1
pca.eigenvalues
=> [1.92592927269225, 0.0740707273077545]
pca.component_matrix
=> GSL::Matrix
[  9.813e-01 
  9.813e-01 ]
pca.communalities
=> [0.962964636346122, 0.962964636346122]

References:

Instance Attribute Summary collapse

Instance Method Summary collapse

Methods included from Summarizable

#summary

Constructor Details

#initialize(matrix, opts = Hash.new) ⇒ PCA

Returns a new instance of PCA.



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
85
86
87
 
# File 'lib/statsample/factor/pca.rb', line 54

def initialize(matrix, opts=Hash.new)
  @use_gsl = opts[:use_gsl]
  opts.delete :use_gsl

  @name=_("Principal Component Analysis")
  @matrix=matrix
  @n_variables=@matrix.column_size      
  @variables_names=(@matrix.respond_to? :fields) ? @matrix.fields : @n_variables.times.map {|i| "VAR_#{i+1}".to_sym }
  
  @matrix_type = @matrix.respond_to?(:_type) ? @matrix._type : :correlation
  
  @m=nil
  
  @rotation_type=Statsample::Factor::Varimax
  
  opts.each{|k,v|
    self.send("#{k}=",v) if self.respond_to? k
  }

  if @use_gsl.nil?
    @use_gsl=Statsample.has_gsl?
  end
  if @matrix.respond_to? :fields
    @variables_names=@matrix.fields
  else
    @variables_names=@n_variables.times.map {|i| "V#{i+1}".to_sym}
  end
  calculate_eigenpairs
  
  if @m.nil?
    # Set number of factors with eigenvalues > 1
    @m=@eigenpairs.find_all {|ev,ec| ev>=1.0}.size
  end
end

Instance Attribute Details

#m ⇒ Object

Number of factors. Set by default to the number of factors with eigen values > 1



44
45
46
 
# File 'lib/statsample/factor/pca.rb', line 44

def m
  @m
end

#matrix_type ⇒ Object

Returns the value of attribute matrix_type.



53
54
55
 
# File 'lib/statsample/factor/pca.rb', line 53

def matrix_type
  @matrix_type
end

#name ⇒ Object

Name of analysis



40
41
42
 
# File 'lib/statsample/factor/pca.rb', line 40

def name
  @name
end

#rotation_type ⇒ Object

Type of rotation. By default, Statsample::Factor::Rotation::Varimax



52
53
54
 
# File 'lib/statsample/factor/pca.rb', line 52

def rotation_type
  @rotation_type
end

#summary_parallel_analysis ⇒ Object

Add to the summary a parallel analysis report



50
51
52
 
# File 'lib/statsample/factor/pca.rb', line 50

def summary_parallel_analysis
  @summary_parallel_analysis
end

#summary_rotation ⇒ Object

Add to the summary a rotation report



48
49
50
 
# File 'lib/statsample/factor/pca.rb', line 48

def summary_rotation
  @summary_rotation
end

#use_gsl ⇒ Object

Use GSL if available



46
47
48
 
# File 'lib/statsample/factor/pca.rb', line 46

def use_gsl
  @use_gsl
end

Instance Method Details

#communalities(m = nil) ⇒ Object 



188
189
190
191
192
193
194
195
196
197
198
199
 
# File 'lib/statsample/factor/pca.rb', line 188

def communalities(m=nil)
  m||=@m
  h=[]
  @n_variables.times do |i|
    sum=0
    m.times do |j|
      sum += (@eigenpairs[j][0].abs*@eigenpairs[j][1][i]**2)
    end
    h.push(sum)
  end
  h
end

#component_matrix(m = nil) ⇒ Object 



145
146
147
148
 
# File 'lib/statsample/factor/pca.rb', line 145

def component_matrix(m=nil)
  var="component_matrix_#{matrix_type}"
  send(var,m)
end

#component_matrix_correlation(m = nil) ⇒ Object

Matrix with correlations between components and variables



170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
 
# File 'lib/statsample/factor/pca.rb', line 170

def component_matrix_correlation(m=nil)
  m||=@m
  raise "m should be > 0" if m<1
  omega_m=::Matrix.build(@n_variables, m) {0}
  gammas=[]
  m.times {|i|
    omega_m.column=i, @eigenpairs[i][1]
    gammas.push(Math::sqrt(@eigenpairs[i][0]))
  }
  gamma_m=::Matrix.diagonal(*gammas)
  cm=(omega_m*(gamma_m)).to_matrix
  
  cm.extend CovariateMatrix
  cm.name=_("Component matrix")
  cm.fields_x = @variables_names
  cm.fields_y = m.times.map { |i| "PC_#{i+1}".to_sym }
  cm
end

#component_matrix_covariance(m = nil) ⇒ Object

Matrix with correlations between components and variables. Based on Härdle & Simar (2003, p.243)



151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
 
# File 'lib/statsample/factor/pca.rb', line 151

def component_matrix_covariance(m=nil)
  m||=@m
  raise "m should be > 0" if m<1
  ff=feature_matrix(m)
  cm=::Matrix.build(@n_variables, m) {0}
  @n_variables.times {|i|
    m.times {|j|
      cm[i,j]=ff[i,j] * Math.sqrt(eigenvalues[j] / @matrix[i,i])
    }
  }
  cm.extend NamedMatrix
  cm.name=_("Component matrix (from covariance)")
  cm.fields_x = @variables_names
  cm.fields_y = m.times.map {|i| "PC_#{i+1}".to_sym }
  
  cm
end

#eigenvalues ⇒ Object

Array with eigenvalues



201
202
203
 
# File 'lib/statsample/factor/pca.rb', line 201

def eigenvalues
  @eigenpairs.collect {|c| c[0] }
end

#eigenvectors ⇒ Object 



204
205
206
207
208
 
# File 'lib/statsample/factor/pca.rb', line 204

def eigenvectors
  @eigenpairs.collect {|c| 
    @use_gsl ? c[1].to_gsl : Daru::Vector.new(c[1])
  }
end

#feature_matrix(m = nil) ⇒ Object

Feature matrix for m factors Returns m eigenvectors as columns. So, i=variable, j=component



106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
 
# File 'lib/statsample/factor/pca.rb', line 106

def feature_matrix(m=nil)
  m||=@m
  if @use_gsl
    omega_m=GSL::Matrix.zeros(@n_variables,m)
    ev=eigenvectors
    m.times do |i|
      omega_m.set_column(i,ev[i])
    end
    omega_m
  else
    omega_m=::Matrix.build(@n_variables, m) {0}
    m.times do |i|
      omega_m.column= i, @eigenpairs[i][1]
    end
    omega_m
  end
end

#principal_components(input, m = nil) ⇒ Object

Returns Principal Components for input matrix or dataset The number of PC to return is equal to parameter m. If m isn’t set, m set to number of PCs selected at object creation. Use covariance matrix



128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
 
# File 'lib/statsample/factor/pca.rb', line 128

def principal_components(input, m=nil)
  if @use_gsl
    data_matrix=input.to_gsl
  else
    data_matrix=input.to_matrix
  end
  m||=@m
  
  raise "data matrix variables<>pca variables" if data_matrix.column_size!=@n_variables
  
  fv=feature_matrix(m)
  pcs=(fv.transpose*data_matrix.transpose).transpose
  
  pcs.extend Statsample::NamedMatrix
  pcs.fields_y = m.times.map { |i| "PC_#{i+1}".to_sym }
  pcs.to_dataframe
end

#report_building(builder) ⇒ Object

:nodoc:



214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
 
# File 'lib/statsample/factor/pca.rb', line 214

def report_building(builder) # :nodoc:
  builder.section(:name=>@name) do |generator|
    generator.text _("Number of factors: %d") % m
    generator.table(:name=>_("Communalities"), :header=>[_("Variable"),_("Initial"),_("Extraction"), _("%")]) do |t|
      communalities(m).each_with_index {|com, i|
        perc=com*100.quo(@matrix[i,i])
        t.row([@variables_names[i], "%0.3f" % @matrix[i,i]  , "%0.3f" % com, "%0.3f" % perc])
      }
    end
    te=total_eigenvalues
    generator.table(:name=>_("Total Variance Explained"), :header=>[_("Component"), _("E.Total"), _("%"), _("Cum. %")]) do |t|
      ac_eigen=0
      eigenvalues.each_with_index {|eigenvalue,i|
        ac_eigen+=eigenvalue
        t.row([_("Component %d") % (i+1), sprintf("%0.3f",eigenvalue), sprintf("%0.3f%%", eigenvalue*100.quo(te)), sprintf("%0.3f",ac_eigen*100.quo(te))])
      }
    end
    
    generator.parse_element(component_matrix(m))
              
    if (summary_rotation)
      generator.parse_element(rotation)
    end
  end
end

#rotation ⇒ Object 



88
89
90
 
# File 'lib/statsample/factor/pca.rb', line 88

def rotation
  @rotation_type.new(component_matrix)
end

#total_eigenvalues ⇒ Object 



91
92
93
 
# File 'lib/statsample/factor/pca.rb', line 91

def total_eigenvalues
  eigenvalues.inject(0) {|ac,v| ac+v}
end