import pandas as pd
import numpy as np

#导入数据
data=pd.read_excel(r'data.xlsx',index_col=0,sheet_name='灰色预测1')
x=np.array(data)
print(x)

#重复代码已省略

#确定参考数列
x_0=np.array([9,9,9,9,8,9,9])

#计算绝对误差
abs_err=abs(x-x_0)
print(abs_err)

#重复代码已省略

#计算两级最小差和两级最大差
min_err_list=[]
max_err_list=[]
for row in abs_err:
    min_err_list.append(min(row))
    max_err_list.append(max(row))
print(min_err_list)
print(max_err_list)
min_err=min(min_err_list)
max_err=max(max_err_list)
print(min_err)
print(max_err)

min_err=min(abs_err.min(axis=1))      #两级最小差
max_err=max(abs_err.max(axis=1))      #两级最大差

#重复代码已省略

#计算关联系数
rho=0.5                         #分辨系数
zeta=np.zeros(abs_err.shape)    #关联系数矩阵

zeta=(self.min_err+self.rho*self.max_err)/(self.abs_err+self.rho*self.max_err)

print(zeta)

#重复代码已省略

#分别计算每个人各个指标关联系数的均值(关联度)
r=[]

for zeta_i in zeta:
    r.append(np.mean(zeta_i))
    
print(r)

import pandas as pd
import numpy as np

#导入数据
data=pd.read_excel(r'data.xlsx',index_col=0,sheet_name='灰色预测1')
x=np.array(data)

#确定参考数列
x_0=np.array([9,9,9,9,8,9,9])

#构建灰色关联类
class GreyRelation:
    #计算关联系数
    def calZeta(self):
        zeta=np.zeros(self.abs_err.shape)        
        zeta=(self.min_err+self.rho*self.max_err)/(self.abs_err+self.rho*self.max_err)
        return zeta
    #计算关联度
    def calR(self):
        r=[]
        for zeta_i in self.zeta:
            r.append(np.mean(zeta_i))
        return r
    def __init__(self,X,X_0,Rho=0.5):
        self.x=X                                        #原始数据
        self.x_0=X_0                                    #参考数据列
        self.abs_err=abs(X-X_0)                         #绝对误差
        self.min_err=min(self.abs_err.min(axis=1))      #两级最小差
        self.max_err=max(self.abs_err.max(axis=1))      #两级最大差
        self.rho=Rho                                    #分辨系数
        self.zeta=self.calZeta()                        #关联系数
        self.r=self.calR()                              #关联度
    def __str__(self):
        s1=f'x:\n{self.x}\n'
        s2=f'x_0:{self.x_0}\n'
        s3=f'abs_err:\n{self.abs_err}\n'
        s4=f'min_err:{self.min_err}\n'
        s5=f'max_err:{self.max_err}\n'
        s6=f'rho:{self.rho}\n'
        s7=f'zeta:\n{self.zeta}\n'
        s8=f'r:{self.r}'
        return(s1+s2+s3+s4+s5+s6+s7+s8)
    
teacher=GreyRelation(x,x_0)
print(teacher)

import numpy as np

#累加生成函数
def cumAddGen(arr,r=1):
    if r==1:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return np.array(res)
    else:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return cumAddGen(np.array(res),r-1)

x_0=np.array([3.2,3.3,3.4,3.6,3.8])

print(cumAddGen(x_0))
print(cumAddGen(cumAddGen(x_0)))
print(cumAddGen(x_0,2))
print(cumAddGen(cumAddGen(cumAddGen(x_0))))
print(cumAddGen(x_0,3))

import numpy as np
from matplotlib import pyplot as plt

#累加生成函数
def cumAddGen(arr,r=1):
    if r==1:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return np.array(res)
    else:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return cumAddGen(np.array(res),r-1)

x_0=np.array([3.2,1.3,4.4,3.6,5.8])
plt.figure(figsize=(8,8),dpi=80)
plt.plot(np.arange(1,6),x_0) 
plt.show()
plt.plot(np.arange(1,6),cumAddGen(x_0)) 
plt.show()

import numpy as np

#累减生成函数
def cumMinusGen(arr,r=1):
    if r==1:
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        return np.array(res)
    else:
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        return cumMinusGen(np.array(res),r-1)
    
x_0=np.array([5,9,14,24,35,46])
print(cumMinusGen(x_0))
print(cumMinusGen(cumMinusGen(x_0)))
print(cumMinusGen(x_0,2))

import numpy as np
import pandas as pd
import math

#导入数据
data=pd.read_excel(r'data.xlsx',index_col=0,sheet_name='灰色预测2')
x_0=np.squeeze(np.array(data).T)    #建立数据时间序列
print(f'x_0={x_0}')

#计算级比
def calClassRatio(arr):
    cr=[]
    for i in range(0,len(arr)-1):
        cr.append(arr[i]/arr[i+1])
    return(np.array(cr))

#判断级比是否落在落在可容覆盖区间内
def checkClassRatio(cr):
    range_pow=2/(len(cr)+2)
    for i in cr:
        if i>=math.pow(math.e,range_pow) or i<=math.pow(math.e,-range_pow):
            return False
    return True

#累加生成函数
def cumAddGen(arr,r=1):
    if r==1:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return np.array(res)
    else:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return cumAddGen(np.array(res),r-1)

print('\n级比检验：')
class_ratio=calClassRatio(x_0)
print(f'λ={class_ratio}')
print(f'λ check res={checkClassRatio(class_ratio)}')

#重复代码已省略

#累加生成函数
def cumAddGen(arr,r=1):
    if r==1:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return np.array(res)
    else:
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        return cumAddGen(np.array(res),r-1)
    
#累减生成函数
def cumMinusGen(arr,r=1):
    if r==1:
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        return np.array(res)
    else:
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        return cumMinusGen(np.array(res),r-1)
    
#构造数据矩阵
def constructDataMatrix(arr_0,arr_1):
    B=np.ones((len(arr_0)-1,2))
    for i in range(0,len(arr_0)-1):
        B[i,0]=-0.5*(arr_1[i]+arr_1[i+1])
    return B

#构造数据向量
def constructDataVector(arr):
    return(np.array([arr[1:]]).T)

#估计发展系数与灰作用量
def estimateU(B,Y):
    return(np.matmul(np.matmul(np.linalg.inv(np.matmul(B.T,B)),B.T),Y))

#对模型灰微分方程进行求解
def solveModel(arr,a,b):
    res=[]
    for k in range(0,len(arr)):
        res.append(((arr[0]-b/a))*math.pow(math.e,-a*k)+b/a)
    return(np.array(res))

print('\nGM(1,1)建模')
x_1=cumAddGen(x_0)              #原始数据一次累加
print(f'x_1={x_1}')

#构造数据矩阵B以及数据向量Y
Y=constructDataVector(x_0)      #数据向量
B=constructDataMatrix(x_0,x_1)  #数据矩阵
print(f'B=\n{B}')
print(f'Y=\n{Y}')

#计算u^
u_estimate=estimateU(B,Y)
a=u_estimate[0,0]
b=u_estimate[1,0]
print(f'u^=\n{u_estimate}')
print(f'a={a}')
print(f'b={b}')

#建立模型并求解
x_1_estimate=solveModel(x_0,a,b)
x_0_estimate=cumMinusGen(x_1_estimate)
print(f'x_1^={x_1_estimate}')
print(f'x_0^={x_0_estimate}')

#重复代码已省略

#计算相对残差
def calRelativeResiduals(arr,arr_e):
    res=[]
    for i in range(0,len(arr)):
        res.append((arr[i]-arr_e[i])/arr[i])
    return(np.array(res))

#计算级比偏差
def calRatioDeviation(cr,a):
    return(1-(1-0.5*a)/(1+0.5*a)*cr)
    
print('\n模型检验：')
relative_residuals=calRelativeResiduals(x_0,x_0_estimate)   #相对残差
ratio_deviation=calRatioDeviation(class_ratio,a)            #级比偏差
print(f'relative_residuals={relative_residuals}')
print(f'ratio_deviation={ratio_deviation}')

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import math

#导入数据
data=pd.read_excel(r'data.xlsx',index_col=0,sheet_name='灰色预测2')
x_0=np.squeeze(np.array(data).T)    #建立数据时间序列

#构造灰色预测类
class GreyModel:
    #级比的计算与检查
    def __calClassRatio(self):
        cr=[]
        for i in range(0,len(self.x_0)-1):
            cr.append(self.x_0[i]/self.x_0[i+1])
        return(np.array(cr))
    def __checkClassRatio(self):
        range_pow=2/(len(self.class_ratio)+2)
        for i in self.class_ratio:
            if i>=math.pow(math.e,range_pow) or i<=math.pow(math.e,-range_pow):
                print("模型级比检查不通过，请检查数据后再试！")
                return False
        return True
    #累加生成函数与累减生成函数
    def __cumAddGen(self,arr,r=1):
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        if r==1:
            return np.array(res)
        else:
            return self.__cumAddGen(arr=np.array(res),r=r-1)   
    def __cumMinusGen(self,arr,r=1):
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        if r==1:            
            return np.array(res)
        else:            
            return self.__cumMinusGen(arr=np.array(res),r=r-1)
    #邻值生成函数
    def __neighborGen(self,arr,a=0.5):
        res=[]
        for i in range(1,len(arr)):
            res.append(a*arr[i]+(1-a)*arr[i-1])
        return np.array(res)
    #构造数据向量Y与数据矩阵B
    def __constructDataVector(self):
        return(np.array([self.x_0[1:]]).T)
    def __constructDataMatrix(self):
        B=np.ones((len(self.x_0)-1,2))
        for i in range(0,len(self.x_0)-1):
            B[i,0]=-self.z_1[i]
        return B
    #估计发展系数与灰作用量
    def __estimateU(self):
        return(np.matmul(np.matmul(np.linalg.inv(np.matmul(self.B.T,self.B)),self.B.T),self.Y))
    #对模型灰微分方程进行求解
    def __solveModel(self):
        res=[]
        for k in range(0,len(self.x_0)):
            res.append(((self.x_0[0]-self.b/self.a))*math.pow(math.e,-self.a*k)+self.b/self.a)
        return(np.array(res))
    #计算相对残差与级比偏差
    def __calRelativeResiduals(self):
        res=[]
        for i in range(0,len(self.x_0)):
            res.append((self.x_0[i]-self.x_0_estimate[i])/self.x_0[i])
        return(np.array(res))
    def __calRatioDeviation(self):
        return(1-(1-0.5*self.a)/(1+0.5*self.a)*self.class_ratio)
    #检查模型
    def __checkModel(self):
        if max(abs(self.relative_residuals))<0.2 and max(abs(self.ratio_deviation))<0.2:
            if max(abs(self.relative_residuals))>=0.1 and max(abs(self.ratio_deviation))>=0.1:
                print("模型精确度仅达到一般要求，认为通过，但还是建议检查数据！")
            return True
        else:
            print("模型精确度不足，请检查数据后重试！")
            return False
    def __init__(self,X_0,weight=0.5):
        self.x_0=X_0                                                    #初始数据列(灰导数)
        self.class_ratio=self.__calClassRatio()                         #级比
        self.cr_accepted=self.__checkClassRatio()                       #级比是否完全落在可容覆盖区间内
        self.x_1=self.__cumAddGen(arr=self.x_0)                         #经过一次累加生成的数据列
        self.weight=weight                                              #权重
        self.z_1=self.__neighborGen(arr=self.x_1,a=self.weight)         #白化背景值
        self.B=self.__constructDataMatrix()                             #数据向量
        self.Y=self.__constructDataVector()                             #数据矩阵
        self.a=self.__estimateU()[0,0]                                  #发展系数
        self.b=self.__estimateU()[1,0]                                  #灰作用量
        self.x_1_estimate=self.__solveModel()                           #x_1的模型估计值
        self.x_0_estimate=self.__cumMinusGen(arr=self.x_1_estimate)     #x_0的模型估计值
        self.relative_residuals=self.__calRelativeResiduals()           #相对残差
        self.ratio_deviation=self.__calRatioDeviation()                 #级比偏差
        self.model_accepted=self.__checkModel()                         #模型精度是否通过检测
    def __str__(self):
        s=''
        s+=f'x_0={self.x_0}\n'
        s+=f'class_ratio={self.class_ratio},accepted={self.cr_accepted}\n'
        s+=f'x_1={self.x_1}\n'
        s+=f'weight={self.weight}\n'
        s+=f'z_1={self.z_1}\n'
        s+=f'B=\n{self.B}\n'
        s+=f'Y=\n{self.Y}\n'
        s+=f'a={self.a},b={self.b}\n'
        s+=f'x_1_estimate={self.x_1_estimate}\n'
        s+=f'x_0_estimate={self.x_0_estimate}\n'
        s+=f'relative_residuals={self.relative_residuals}\n'
        s+=f'ratio_deviation={self.ratio_deviation}\n'
        s+=f'model_accepted={self.model_accepted}'
        return s
    #接口：通过图像检查原始值与模型预测值之间的关系
    def showModel(self):
        #原始值图像曲线
        x1=np.arange(1,len(self.x_0)+1)
        y1=self.x_0
        #模型预测值图像曲线
        x2=np.arange(1,len(self.x_0_estimate)+1)
        y2=self.x_0_estimate
        plt.plot(x1,y1,x2,y2,marker='o')                       
        plt.xlabel("K - Time")
        plt.ylabel("Value")
        plt.legend(labels=["Estimated Value","Real Value"],loc="lower right") 
        plt.show()
    #接口：通过模型进行预测,返回从k=1开始直到k=target的完整预测结果数列
    def forecast(self,target):
        res=[]
        for k in range(0,target):
            res.append(((self.x_0[0]-self.b/self.a))*math.pow(math.e,-self.a*k)+self.b/self.a)
        return(self.__cumMinusGen(arr=(np.array(res))))
    #接口：通过模型进行预测,并绘制预测结果图
    def forecastAndShow(self,start=1,end=None):
        forecast_data=self.forecast(target=end)
        x=np.arange(start,end+1)
        y=forecast_data
        plt.plot(x,y,marker='o')
        plt.xlabel("K - Time")
        plt.ylabel("Value")
        for a,b in zip(x,y):
            plt.text(a,b+0.02,'%.4f'%b,ha='center',va='bottom')
        plt.show()
            
A=GreyModel(x_0)
print(A)
A.showModel()
print(A.forecast(target=12))
A.forecastAndShow(end=12)

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
import math

file='./data/灰色预测_SARS.xlsx'

retail_data_read=pd.read_excel(file,sheet_name='商品零售额',index_col=0)
tourist_data_read=pd.read_excel(file,sheet_name='接待海外旅游人数',index_col=0)
service_data_read=pd.read_excel(file,sheet_name='综合服务业累计数据',index_col=0)

print(f'表1：商品的零售额（单位：亿元）\n{retail_data_read}\n\n')
print(f'表2 接待海外旅游人数（单位：万人）\n{tourist_data_read}\n\n')
print(f'表3 综合服务业累计数据（单位：亿元）\n{service_data_read}\n\n')

#剔除了2003年的数据
retail_data=np.array(retail_data_read[:-1])      #零售业数据
tourist_data=np.array(tourist_data_read[:-1])    #旅游业数据
service_data=np.array(service_data_read[:-1])    #服务业数据

#2003年的数据
retail_2003=np.squeeze(np.array(retail_data_read[-1:]))
tourist_2003=np.squeeze(np.array(tourist_data_read[-1:]))
service_2003=np.squeeze(np.array(service_data_read[-1:]))

print(tourist_data)
print(tourist_data)
print(service_data)
print(tourist_2003)
print(tourist_2003)
print(service_2003)

#构造灰色预测类
class GreyModel:
    #级比的计算与检查
    def __calClassRatio(self):
        cr=[]
        for i in range(0,len(self.x_0)-1):
            cr.append(self.x_0[i]/self.x_0[i+1])
        return(np.array(cr))
    def __checkClassRatio(self):
        range_pow=2/(len(self.class_ratio)+2)
        for i in self.class_ratio:
            if i>=math.pow(math.e,range_pow) or i<=math.pow(math.e,-range_pow):
                print("模型级比检查不通过，请检查数据后再试！")
                return False
        return True
    #累加生成函数与累减生成函数
    def __cumAddGen(self,arr,r=1):
        res=[]
        for i in range(1,len(arr)+1):
            res.append(sum(arr[0:i]))
        if r==1:
            return np.array(res)
        else:
            return self.__cumAddGen(arr=np.array(res),r=r-1)   
    def __cumMinusGen(self,arr,r=1):
        res=[arr[0]]
        for i in range(1,len(arr)):
            res.append(arr[i]-arr[i-1])
        if r==1:            
            return np.array(res)
        else:            
            return self.__cumMinusGen(arr=np.array(res),r=r-1)
    #邻值生成函数
    def __neighborGen(self,arr,a=0.5):
        res=[]
        for i in range(1,len(arr)):
            res.append(a*arr[i]+(1-a)*arr[i-1])
        return np.array(res)
    #构造数据向量Y与数据矩阵B
    def __constructDataVector(self):
        return(np.array([self.x_0[1:]]).T)
    def __constructDataMatrix(self):
        B=np.ones((len(self.x_0)-1,2))
        for i in range(0,len(self.x_0)-1):
            B[i,0]=-self.z_1[i]
        return B
    #估计发展系数与灰作用量
    def __estimateU(self):
        return(np.matmul(np.matmul(np.linalg.inv(np.matmul(self.B.T,self.B)),self.B.T),self.Y))
    #对模型灰微分方程进行求解
    def __solveModel(self):
        res=[]
        for k in range(0,len(self.x_0)):
            res.append(((self.x_0[0]-self.b/self.a))*math.pow(math.e,-self.a*k)+self.b/self.a)
        return(np.array(res))
    #计算相对残差与级比偏差
    def __calRelativeResiduals(self):
        res=[]
        for i in range(0,len(self.x_0)):
            res.append((self.x_0[i]-self.x_0_estimate[i])/self.x_0[i])
        return(np.array(res))
    def __calRatioDeviation(self):
        return(1-(1-0.5*self.a)/(1+0.5*self.a)*self.class_ratio)
    #检查模型
    def __checkModel(self):
        if max(abs(self.relative_residuals))<0.2 and max(abs(self.ratio_deviation))<0.2:
            if max(abs(self.relative_residuals))>=0.1 and max(abs(self.ratio_deviation))>=0.1:
                print("模型精确度仅达到一般要求，认为通过，但还是建议检查数据！")
            return True
        else:
            print("模型精确度不足，请检查数据后重试！")
            return False
    def __init__(self,X_0,weight=0.5):
        self.x_0=X_0                                                    #初始数据列(灰导数)
        self.class_ratio=self.__calClassRatio()                         #级比
        self.cr_accepted=self.__checkClassRatio()                       #级比是否完全落在可容覆盖区间内
        self.x_1=self.__cumAddGen(arr=self.x_0)                         #经过一次累加生成的数据列
        self.weight=weight                                              #权重
        self.z_1=self.__neighborGen(arr=self.x_1,a=self.weight)         #白化背景值
        self.B=self.__constructDataMatrix()                             #数据向量
        self.Y=self.__constructDataVector()                             #数据矩阵
        self.a=self.__estimateU()[0,0]                                  #发展系数
        self.b=self.__estimateU()[1,0]                                  #灰作用量
        self.x_1_estimate=self.__solveModel()                           #x_1的模型估计值
        self.x_0_estimate=self.__cumMinusGen(arr=self.x_1_estimate)     #x_0的模型估计值
        self.relative_residuals=self.__calRelativeResiduals()           #相对残差
        self.ratio_deviation=self.__calRatioDeviation()                 #级比偏差
        self.model_accepted=self.__checkModel()                         #模型精度是否通过检测
    def __str__(self):
        s=''
        s+=f'x_0={self.x_0}\n'
        s+=f'class_ratio={self.class_ratio},accepted={self.cr_accepted}\n'
        s+=f'x_1={self.x_1}\n'
        s+=f'weight={self.weight}\n'
        s+=f'z_1={self.z_1}\n'
        s+=f'B=\n{self.B}\n'
        s+=f'Y=\n{self.Y}\n'
        s+=f'a={self.a},b={self.b}\n'
        s+=f'x_1_estimate={self.x_1_estimate}\n'
        s+=f'x_0_estimate={self.x_0_estimate}\n'
        s+=f'relative_residuals={self.relative_residuals}\n'
        s+=f'ratio_deviation={self.ratio_deviation}\n'
        s+=f'model_accepted={self.model_accepted}'
        return s
    #接口：通过图像检查原始值与模型预测值之间的关系
    def showModel(self):
        #原始值图像曲线
        x1=np.arange(1,len(self.x_0)+1)
        y1=self.x_0
        #模型预测值图像曲线
        x2=np.arange(1,len(self.x_0_estimate)+1)
        y2=self.x_0_estimate
        plt.plot(x1,y1,x2,y2,marker='o')                       
        plt.xlabel("K - Time")
        plt.ylabel("Value")
        plt.legend(labels=["Estimated Value","Real Value"],loc="lower right") 
        plt.show()
    #接口：通过模型进行预测,返回从k=1开始直到k=target的完整预测结果数列
    def forecast(self,target):
        res=[]
        for k in range(0,target):
            res.append(((self.x_0[0]-self.b/self.a))*math.pow(math.e,-self.a*k)+self.b/self.a)
        return(self.__cumMinusGen(arr=(np.array(res))))
    #接口：通过模型进行预测,并绘制预测结果图
    def forecastAndShow(self,start=1,end=None):
        forecast_data=self.forecast(target=end)
        x=np.arange(start,end+1)
        y=forecast_data
        plt.plot(x,y,marker='o')
        plt.xlabel("K - Time")
        plt.ylabel("Value")
        for a,b in zip(x,y):
            plt.text(a,b+0.02,'%.4f'%b,ha='center',va='bottom')
        plt.show()

#计算1997-2002年各行业的年均值
retail_year_avg=np.average(retail_data,axis=1)          #零售业年均值
tourist_year_avg=np.average(tourist_data,axis=1)        #旅游业年均值
service_year_avg=np.average(service_data,axis=1)        #服务业年均值

#计算各行业每月在年内的总值占比
retail_month_sum=np.sum(retail_data,axis=0)             #零售业月总值
tourist_month_sum=np.sum(tourist_data,axis=0)           #旅游业月总值
service_month_sum=np.sum(service_data,axis=0)           #服务业月总值
retail_share=retail_month_sum/np.sum(retail_data)       #零售业月占比
tourist_share=tourist_month_sum/np.sum(tourist_data)    #旅游业月占比
service_share=service_month_sum/np.sum(service_data)    #服务业月占比

print('各行业年均值：')
print(retail_year_avg)
print(tourist_year_avg)
print(service_year_avg)
print('各行业月占比：')
print(retail_share)
print(tourist_share)
print(service_share)

retail=GreyModel(retail_year_avg,weight=0.4)   #零售业模型,权重根据提示取0.4
print(retail)
retail.showModel()

retail_2003_est_sum=retail.forecast(target=len(retail.x_0)+1)[len(retail.x_0)]*12   #2003年预计总值
retail_2003_est=retail_2003_est_sum*retail_share                                    #2003年预计月分布值

print(f'2003年预计零售业总值为：{retail_2003_est_sum}')
print('按月分布如下：')
print(retail_2003_est)

#时间横轴
x=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEPT','OCT','NOV','DEC']

#数据纵轴
y1=retail_2003_est              #预测值
y2=retail_2003                  #实际值

plt.figure(figsize=(16,8),dpi=80)
plt.title("Retail Compared by Month")
plt.xlabel("Month")
plt.ylabel("Retail")
plt.plot(x,y1,x,y2,marker='o')
plt.legend(labels=["Estimated Value","Real Value"],loc="lower right")  
for a,b,c in zip(x,y1,y2):
    plt.text(a,b,'%.1f'%b,ha='center',va='bottom')
    plt.text(a,c,'%.1f'%c,ha='center',va='bottom')    
plt.show()

tourist=GreyModel(tourist_year_avg,weight=0.5)   #旅游业模型,权重根据提示取0.5
print(tourist)
tourist.showModel()

tourist_2003_est_sum=tourist.forecast(target=len(tourist.x_0)+1)[len(tourist.x_0)]*12   #2003年预计总值
tourist_2003_est=tourist_2003_est_sum*tourist_share                                     #2003年预计月分布值

print(f'2003年预计旅游业总值为：{tourist_2003_est_sum}')
print('按月分布如下：')
print(tourist_2003_est)

#时间横轴
x=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEPT','OCT','NOV','DEC']

#数据纵轴
y1=tourist_2003_est              #预测值
y2=tourist_2003                  #实际值

plt.figure(figsize=(16,8),dpi=80)
plt.title("Tourist Compared by Month")
plt.xlabel("Month")
plt.ylabel("Tourist")
plt.plot(x,y1,x,y2,marker='o')
plt.legend(labels=["Estimated Value","Real Value"],loc="lower right")  
for a,b,c in zip(x,y1,y2):
    plt.text(a,b,'%.1f'%b,ha='center',va='bottom')
    plt.text(a,c,'%.1f'%c,ha='center',va='bottom')    
plt.show()

service=GreyModel(service_year_avg,weight=0.5)   #旅游业模型,权重根据提示取0.5
print(service)
service.showModel()

service_2003_est_sum=service.forecast(target=len(service.x_0)+1)[len(service.x_0)]      #2003年预计均值
service_2003_est=12*service_2003_est_sum*service_share                                  #2003年预计月分布值

print(f'2003年预计服务业总值为：{service_2003_est_sum}')
print('按月分布如下：')
print(service_2003_est)

#时间横轴
x=['FEB','MAR','APR','MAY','JUN','JUL','AUG','SEPT','OCT','NOV','DEC']

#数据纵轴
y1=service_2003_est              #预测值
y2=service_2003                  #实际值

plt.figure(figsize=(16,8),dpi=80)
plt.title("Service Compared by Month")
plt.xlabel("Month")
plt.ylabel("Service")
plt.plot(x,y1,x,y2,marker='o')
plt.legend(labels=["Estimated Value","Real Value"],loc="lower right")  
for a,b,c in zip(x,y1,y2):
    plt.text(a,b,'%.1f'%b,ha='center',va='bottom')
    plt.text(a,c,'%.1f'%c,ha='center',va='bottom')    
plt.show()