# -*- coding: UTF-8 -*-

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from statsmodels.stats.api import anova_lm
import matplotlib.pyplot as plt
import pandas as pd
from patsy import dmatrices
import itertools as it
import random


# Load data 讀取數(shù)據(jù)
df = pd.read_csv('data3.1.csv',encoding='gbk')
print(df)


target = 'y'
variate = set(df.columns) #獲取列名
variate.remove(target) #去除無關(guān)列
variate.remove('地區(qū)')

#定義多個(gè)數(shù)組，用來分別用來添加變量，刪除變量
x = []
variate_add = []
variate_del = variate.copy()
# print(variate_del)
y = random.sample(variate,3) #隨機(jī)生成一個(gè)選模型，3為變量的個(gè)數(shù)
print(y)
#將隨機(jī)生成的三個(gè)變量分別輸入到 添加變量和刪除變量的數(shù)組
for i in y:
 variate_add.append(i)
 x.append(i)
 variate_del.remove(i)

global aic #設(shè)置全局變量 這里選擇AIC值作為指標(biāo)
formula="{}~{}".format("y","+".join(variate_add)) #將自變量名連接起來
aic=smf.ols(formula=formula,data=df).fit().aic #獲取隨機(jī)函數(shù)的AIC值，與后面的進(jìn)行對比
print("隨機(jī)化選模型為：{}~{}，對應(yīng)的AIC值為：{}".format("y","+".join(variate_add), aic))
print("\n")



#添加變量
def forwark():
 score_add = []
 global best_add_score
 global best_add_c
 print("添加變量")
 for c in variate_del:
  formula = "{}~{}".format("y", "+".join(variate_add+[c]))
  score = smf.ols(formula = formula, data = df).fit().aic
  score_add.append((score, c)) #將添加的變量，以及新的AIC值一起存儲在數(shù)組中
  
  print('自變量為{}，對應(yīng)的AIC值為：{}'.format("+".join(variate_add+[c]), score))

 score_add.sort(reverse=True) #對數(shù)組內(nèi)的數(shù)據(jù)進(jìn)行排序，選擇出AIC值最小的
 best_add_score, best_add_c = score_add.pop()
 
 print("最小AIC值為：{}".format(best_add_score))
 print("\n")

#刪除變量
def back():
 score_del = []
 global best_del_score
 global best_del_c
 print("剔除變量")
 for i in x:

  select = x.copy() #copy一個(gè)集合，避免重復(fù)修改到原集合
  select.remove(i)
  formula = "{}~{}".format("y","+".join(select))
  score = smf.ols(formula = formula, data = df).fit().aic
  print('自變量為{}，對應(yīng)的AIC值為：{}'.format("+".join(select), score))
  score_del.append((score, i))

 score_del.sort(reverse=True) #排序，方便將最小值輸出
 best_del_score, best_del_c = score_del.pop() #將最小的AIC值以及對應(yīng)剔除的變量分別賦值
 print("最小AIC值為：{}".format(best_del_score))
 print("\n")

print("剩余變量為：{}".format(variate_del))
forwark()
back()

while variate:
  
#  forwark()
#  back()
 if(aic < best_add_score < best_del_score or aic < best_del_score < best_add_score):
  print("當(dāng)前回歸方程為最優(yōu)回歸方程，為{}~{}，AIC值為：{}".format("y","+".join(variate_add), aic))
  break
 elif(best_add_score < best_del_score < aic or best_add_score < aic < best_del_score):
  print("目前最小的aic值為{}".format(best_add_score))
  print('選擇自變量：{}'.format("+".join(variate_add + [best_add_c]))) 
  print('\n')
  variate_del.remove(best_add_c)
  variate_add.append(best_add_c)
  print("剩余變量為：{}".format(variate_del))
  aic = best_add_score
  forwark()
 else:
  print('當(dāng)前最小AIC值為：{}'.format(best_del_score))
  print('需要剔除的變量為：{}'.format(best_del_c))
  aic = best_del_score #將AIC值較小的選模型AIC值賦給aic再接著下一輪的對比
  x.remove(best_del_c) #在原集合上剔除選模型所對應(yīng)剔除的變量
  back()