import numpy as np
import pandas as pd
#数据读取
#利用pandas读取csv，读取的数据为DataFrame对象
data = pd.read_csv('sn.csv')
# 将DataFrame对象转化为数组,数组的最后一列为预报对象
data= data.values.copy()
# print(data)
 
# 计算回归系数，参数
def get_regre_coef(X,Y):
  S_xy=0
  S_xx=0
  S_yy=0
  # 计算预报因子和预报对象的均值
  X_mean = np.mean(X)
  Y_mean = np.mean(Y)
  for i in range(len(X)):
    S_xy += (X[i] - X_mean) * (Y[i] - Y_mean)
    S_xx += pow(X[i] - X_mean, 2)
    S_yy += pow(Y[i] - Y_mean, 2)
  return S_xy/pow(S_xx*S_yy,0.5)
#构建原始增广矩阵
def get_original_matrix():
  # 创建一个数组存储相关系数,data.shape几行(维)几列，结果用一个tuple表示
  # print(data.shape[1])
  col=data.shape[1]
  # print(col)
  r=np.ones((col,col))#np.ones参数为一个元组(tuple)
  # print(np.ones((col,col)))
  # for row in data.T:#运用数组的迭代，只能迭代行，迭代转置后的数组，结果再进行转置就相当于迭代了每一列
    # print(row.T)
  for i in range(col):
    for j in range(col):
      r[i,j]=get_regre_coef(data[:,i],data[:,j])
  return r

def get_vari_contri(r):
  col = data.shape[1]
   #创建一个矩阵来存储方差贡献值
  v=np.ones((1,col-1))
  # print(v)
  for i in range(col-1):
    # v[0,i]=pow(r[i,col-1],2)/r[i,i]
    v[0, i] = pow(r[i, col - 1], 2) / r[i, i]
  return v
#选择因子是否进入方程，
#参数说明：r为增广矩阵，v为方差贡献值，k为方差贡献值最大的因子下标,p为当前进入方程的因子数
def select_factor(r,v,k,p):
  row=data.shape[0]#样本容量
  col=data.shape[1]-1#预报因子数
  #计算方差比
  f=(row-p-2)*v[0,k-1]/(r[col,col]-v[0,k-1])
  # print(calc_vari_contri(r))
  return f

#计算第零步增广矩阵
r=get_original_matrix()
# print(r)
#计算方差贡献值
v=get_vari_contri(r)
print(v)
#计算方差比

# #计算方差比
# print(data.shape[0])
f=select_factor(r,v,4,0)
print(f)
#######输出##########
22.79852020138227

#逐步回归分析与计算
#通过矩阵转换公式来计算各部分增广矩阵的元素值
def convert_matrix(r,k):
  col=data.shape[1]
  k=k-1#从第零行开始计数
  #第k行的元素单不属于k列的元素
  r1 = np.ones((col, col)) # np.ones参数为一个元组(tuple)
  for i in range(col):
    for j in range(col):
      if (i==k and j!=k):
        r1[i,j]=r[k,j]/r[k,k]
      elif (i!=k and j!=k):
        r1[i,j]=r[i,j]-r[i,k]*r[k,j]/r[k,k]
      elif (i!= k and j== k):
        r1[i,j] = -r[i,k]/r[k,k]
      else:
        r1[i,j] = 1/r[k,k]
  return r1

#矩阵转换，计算第一步矩阵
r=convert_matrix(r,4)
# print(r)
#计算第一步方差贡献值
v=get_vari_contri(r)
#print(v)
f=select_factor(r,v,1,1)
print(f)
#########输出#####
108.22390933074443

f=select_factor(r,v,2,2)
print(f)
##########输出#########
5.025864648951804

#选择是否剔除因子，
#参数说明：r为增广矩阵，v为方差贡献值，k为方差贡献值最大的因子下标,t为当前进入方程的因子数
def delete_factor(r,v,k,t):
  row = data.shape[0] # 样本容量
  col = data.shape[1] - 1 # 预报因子数
  # 计算方差比
  f = (row - t - 1) * v[0, k - 1] / r[col, col]
  # print(calc_vari_contri(r))
  return f
#因子3的引进检验值0.018233473487350636
f=select_factor(r,v,3,3)
print(f)
#因子4的剔除检验值1.863262422188088
f=delete_factor(r,v,4,3)
print(f)

#因子4的引进检验值1.8632624221880876，小于3.28不能引进
f=select_factor(r,v,4,2)
print(f)
#因子1的剔除检验值146.52265486251397,大于3.28不能剔除
f=delete_factor(r,v,1,2)
print(f)

#因子1和因子2进入方程
#对进入方程的预报因子进行多元回归
# regs=LinearRegression()
X=data[:,0:2]
Y=data[:,4]
X=np.mat(np.c_[np.ones(X.shape[0]),X])#为系数矩阵增加常数项系数
Y=np.mat(Y)#数组转化为矩阵
#print(X)
B=np.linalg.inv(X.T*X)*(X.T)*(Y.T)
print(B.T)#输出系数,第一项为常数项，其他为回归系数
###输出##
#[[52.57734888 1.46830574 0.66225049]]