# 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