编辑代码

# AUthor:Housyou 学院:地空学院

LOGICAL_0 = 1e-10  # 逻辑0值
THRESHOLD = 0.85  # 主成分分析的累计百分比阈值

class Array  # 为内置的数组Array类添加几个方法
  def *(array)  # 数组是个矩阵时的矩阵乘法
    res = nil  # 将要返回的结果，新的矩阵
    if self[0].size == array.size && !array.empty?  # 如果矩阵形状可以进行乘法
      dim = array[0].size  # 后矩阵的列数
      res = Array.new(self.size){Array.new(dim)}
      self.each_with_index do |line,i|  #对前矩阵每行进行遍历，line是每一行，i是索引
        dim.times do |j|  # 对后矩阵每一列进行遍历，j是索引
          value = 0  # 行乘列的结果
          line.each_with_index do |a,k|  # a是前矩阵一行中的每个元素，k是索引
            value += a * array[k][j]  # 行乘列累计求和
          end
          res[i][j] = value  # 将结果保存
        end
      end
    end
    res  # 将矩阵乘法的结果返回，Ruby会自动将最后一次执行的结果作为返回值
  end

  def trans  # 转置
    res = Array.new(self[0].size){Array.new(self.size)}  # 结果的行数和列数对换
    self.each_with_index do |line,i|  # 遍历赋值
      line.each_with_index do |num,j|
        res[j][i] = num
      end
    end
    res  # 返回结果
  end

  def sort_zip(array)  # 对自身进行排序，顺带将另一个数组array按自己的顺序排序
    a1 = self.clone  # 拷贝副本，防止改变自身结果
    a2 = array.clone
    for i in 1...a1.size  # 直接插入排序，因为规模比较小所以不需要用高效的排序算法
      key = a1[i]
      ano = a2[i]
      j = i - 1
      while j >= 0 && yield(key,a1[j])  # yield是调用方法时传入的排序原则
        a1[j + 1] = a1[j]
        a2[j + 1] = a2[j]
        j -= 1
      end
      a1[j + 1] = key
      a2[j + 1] = ano
    end
    {:self=>a1,:other=>a2}  # 排序的结果用散列记录，:self键返回自身排序结果，另一个array排序的结果的键是:other
  end

  def formal  # 矩阵转化为标准化的字符串形式
    res = ""
    self.each do |line|  # 遍历每行
      res += "|"
      line.each do |num|  # 遍历每列
        res += num.round(3).to_s + "\t"
      end
      res += "|\n"
    end
    res
  end

  def find_max  # 找到绝对值最大的元素，返回值和索引
    max_num = 0  # 记录最大的值
    pos = nil
    self.each_with_index do |line,i|
      line.each_with_index do |num,j|
        if i != j
          abs = num.abs
          if abs > max_num
            pos = [i,j]
            max_num = abs
          end
        end
      end
    end
    {value:max_num,i:pos[0],j:pos[1]}
  end
end

def standarlize(address)  # 输出标准化的数据

  def get_ori_data(address)  # 读取原始数据
    res = Array.new  # 保存数据的矩阵数组
    File.open(address) do |file|  # 读取文件
      file.each_line do |line|  # 按行读取
        res.push(line.chomp.split("\t").map {|i| i.to_f})  # 切割后转换为数组，保存数组
      end
    end
    res  # 返回结果
  end

  def standarlize_a_line(array)  # 对一组x进行标准化
    sum,sum2 = 0.0,0.0  # 累计和、累计平方和
    array.each do |num|  # 计算累计和和累计平方和
      sum += num
      sum2 += num * num
    end
    mean = sum / array.size  # 平均数
    s = Math.sqrt((sum2 - sum * mean) / (array.size - 1))  # 标准差
    if s > 1e-10  # 标准差不为0，就输出标准化的数据
      array.map do |num|  # map函数直接对每个数进行替换
        ((num - mean) / s)
      end
    else  # 如果标准差为0
      array.map do |num|
        0  # 平均数也是0，所以每个数都是0
      end
    end
  end

  (get_ori_data(address).trans).map{|line| standarlize_a_line(line)}  # 对整个原始数据进行标准化，返回结果
end

def r_array(data)  # 从标准化数据计算相关系数矩阵
  dim = data.size
  free_degree = data[0].size - 1
  res = Array.new(dim){Array.new(dim)}  # 记录结果
  dim.times do |i|  # 遍历每组x
    dim.times do |j|  # 再次遍历每组x
      r = 1  # 如果 i 和 j 相等，相关系数就是1
      if i > j  # 如果是下半矩阵，直接复制对称位置的值
        r = res[j][i]
      elsif i < j  # 如果是上半矩阵就计算相关系数
        r = 0
        data[i].zip(data[j]) do |a,b|  # 遍历同一个索引对应的两组x的数
          r += a * b  # 计算乘积和
        end
        r /= free_degree  # 相关系数还要除以自由度，因为之前计算方差时是用自由度计算的
      end
      res[i][j] = r  # 将相关系数赋到对应的位置
    end
  end
  res
end

def jacobi(r)  # 雅可比法计算特征向量和特征值
  def get_i(n)  # 获得一个 n 阶单位矩阵
    array = Array.new(n){Array.new(n){0}}  # 对角线外都是0
    n.times {|k| array[k][k] = 1}  # 对角线上是1
    array  # 返回 array
  end
  dim = r.size  # 是方阵
  t = get_i(dim)
  loop do  # 不算到非对角元素全为0就不停止循环
    res = r.find_max  # 找到最大的元素及其索引
    if res[:value] > LOGICAL_0  # 如果非对角元素有不为0的
      i,j = res[:i],res[:j]
      delta = r[i][i] - r[j][j]
      theta = delta.abs < LOGICAL_0 ? -Math::PI/4 : Math.atan(2*r[i][j]/delta)/2  # 取 atan ∞ = -90°
      new_t = get_i(dim)  # 本次迭代的变换矩阵
      cos = Math.cos(theta)
      sin = Math.sin(theta)
      new_t[i][i],new_t[j][j] = cos,cos
      new_t[i][j],new_t[j][i] = -sin,sin
      r = new_t.trans*r*new_t
      r[i][j] = 0  # 为了提高精确度，直接令消去的位置为0
      t *= new_t  # 连乘最后争取获得特征向量
    else  # 满足条件后结束循环
      values = Array.new(dim)  # 记录特征值
      dim.times do |i|
        values[i] = r[i][i]
      end
      return values.sort_zip(t.trans){|a,b| a > b}  # 按特征值对特征值和特征向量排序后返回
    end
  end
end

def contribution(values)  # 计算贡献率
  n = values.size
  res = Array.new(n)
  sum = 0
  values.each {|v| sum += v}  # 计算特征值之和
  n.times {|i| res[i] = values[i] / sum}  # 计算贡献率
  res  # 返回结果
end

def sum_contribution(con)  # 计算累计贡献率及累计贡献率超过阈值时的主成分数量
  thres_index = 0  # 阈值对应的主成分索引
  res = Array.new(con.size){0}  # 计算结果初始化为0
  con.each_with_index do |rate,index|  # 遍历特征值
    res[index] = res[index - 1] + rate  # 如果是第一个特征值，会在最后一个特征值的累计贡献率（目前是0）基础上相加
    if thres_index == 0 && res[index] > THRESHOLD # 超过阈值时记录索引
      thres_index = index
    end
  end
  {:threshold=>thres_index,sum_con:res}  # 返回索引和累计贡献率
end

def load_array(vectors,values)  # 计算载荷矩阵
  dim = [vectors.size,vectors[0].size]
  res = Array.new(dim[1]){Array.new(dim[0])}
  dim[1].times do |i|
    dim[0].times do |j|
      res[i][j] = Math.sqrt(values[j]) * vectors[j][i]
    end
  end
  res
end

def scores(vectors,data)  # 计算主成分得分矩阵
  (vectors * data).trans
end

def analysis(scores_array,num_of_component=1,amount=4,highest=true)  # 分析主成分得分，返回流域索引（从1开始）
  # scores_array是得分矩阵，num_of_component是第几主成分，amount是最高或最低的选取数量，highest判断是否选取最高的几个
  nc = num_of_component - 1
  dim = scores_array.size
  order = Array.new(dim)  # 排序前的索引
  dim.times {|i| order[i] = i + 1}  # 就是原顺序，再加上从1开始计数
  order = (scores_array.sort_zip(order){|a,b| a[nc] > b[nc]})[:other]  # 按主成分得分排序后取order的排序结果
  res = Array.new(amount)  # 待返回的流域索引（从1开始）
  if highest  # 若想最高的，就从头添加
    amount.times do |i|
      res[i] = order[i]
    end
  else  # 反之则从尾部开始添加
    amount.times do |i|
      res[i] = order[-i-1]
    end
  end
  res
end

def main
  data = standarlize("input.txt")  # 获得标准化处理后的数据
  cca = r_array(data)  # 相关系数矩阵

  puts "相关系数矩阵"
  puts cca.formal

  jac = jacobi(cca)
  con = contribution(jac[:self])
  sc = sum_contribution(con)

  puts "特征值及其百分率"
  puts "主成分\t特征值\t百分数\t累计百分数"
  jac[:self].each_with_index do |v,i|
    puts  "#{i+1}\t#{v.round(3)}\t#{con[i].round(3)}\t#{sc[:sum_con][i].round(4)}"
  end

  vectors = jac[:other][0..sc[:threshold]]
  puts "主成分载荷矩阵"
  puts load_array(vectors,jac[:self]).formal

  s = scores(jac[:other],data)
  puts "第一主成分得分最高的流域"
  res = analysis(s)
  res.each do |i|
    print i.to_s + "\t"
  end
  print "\n"
  res.each do |i|
    print s[i - 1][0].round(3).to_s + "\t"
  end
  print "\n"

  puts "第一主成分得分最低的流域"
  res = analysis(s,1,4,false)
  res.each do |i|
    print i.to_s + "\t"
  end
  print "\n"
  res.each do |i|
    print s[i - 1][0].round(3).to_s + "\t"
  end
end

main