class TimingModel:
    def getTimeSeries(self,ts):
        if type(ts)!=type(None):
            return ts
        else:
            return np.arange(0,self.t+1)
    def __init__(self,arr:np.ndarray,time_series=None) -> None:
        self.arr=arr/1.                                     # 时间序列数据
        self.t=self.arr.size                                # 时间序列长度
        self.time_series=self.getTimeSeries(time_series)    # 时间序列

class SimpleMovingAverage(TimingModel):    
    def simpleMovingAverage(self):
        M=np.array([],dtype=np.float64)
        t=0
        while t<self.t+1:
            M=np.insert(M,M.size,np.sum(self.arr[t-self.N:t])/self.N)
            t+=1
        M[M==0]=np.nan
        return M
    def calStdError(self):
        return np.sqrt(np.sum(np.power(self.M[self.N:-1]-self.arr[self.N:],2))/(self.t-self.N))
    def __init__(self,arr:np.ndarray,N,time_series=None) -> None:
        super().__init__(arr,time_series)
        self.N=N                                            # 移动平均的项数
        self.M=self.simpleMovingAverage()                   # 一次预测数列
        self.Mt=self.M[-1:][0]                              # 一次预测值
        self.std_err=self.calStdError()                     # 预测标准误差                           
    def show(self,lite=False):
        print(f'N={self.N}的预测结果如下：')
        raw_data=np.insert(self.arr,self.t,np.nan).reshape(1,self.t+1)
        est_data=self.M.reshape(1,self.M.size)
        columns=self.time_series
        index=['观测值','预测值']
        data=np.r_[raw_data,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'预测结果为：{round(self.Mt,4)}，标准误差为{round(self.std_err,4)}')

data=np.array([533.8,574.6,606.9,649.8,705.1,722.0,816.4,892.7,963.9,1015.1,1102.7])
ts=['JAN','FEB','MAR','APR','MAY','JUN','JUL','AUG','SEPT','OCT','NOV','DEC']
profit=SimpleMovingAverage(data,4,time_series=ts)
profit.show()

class WeightedMovingAverage(TimingModel):
    def weightedMovingAverage(self):
        M=np.zeros(self.N)
        t=self.N
        while t<self.t+1:
            M=np.insert(M,M.size,(np.sum(np.flip(self.arr[t-self.N:t])*self.weight))/np.sum(self.weight))
            t+=1
        M[M==0]=np.nan
        return M
    def calRelativeError(self):
        return (self.arr-self.M[:-1])/self.arr
    def __init__(self,arr:np.ndarray,weight:np.ndarray,time_series=None) -> None:
        super().__init__(arr,time_series)
        self.weight=weight                                  # 分配权重
        self.N=self.weight.size                             # 权重长度
        self.M=self.weightedMovingAverage()                 # 一次预测数列
        self.Mt=self.M[-1:][0]                              # 一次预测结果
        self.relative_err=self.calRelativeError()           # 相对误差
        self.avg_rel_err=np.nanmean(self.relative_err)      # 平均相对误差
        self.adj_Mt=self.Mt/(1-self.avg_rel_err)            # 修正预测结果
    def show(self,lite=False):
        print(f'取weight={self.weight}的预测结果：')
        raw_data=np.insert(self.arr,self.t,np.nan).reshape(1,self.t+1)
        est_data=self.M.reshape(1,self.M.size)
        rel_err=np.insert(self.relative_err,self.t,np.nan).reshape(1,self.t+1)
        data=np.r_[raw_data,est_data,rel_err*100]
        columns=self.time_series
        index=['观测值','预测值','相对误差%']
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'平均相对误差={round(self.avg_rel_err*100,4)}%,修正预测值={round(self.adj_Mt,4)}')

data=np.array([6.35,6.20,6.22,6.66,7.15,7.89,8.72,8.94,9.28,9.8])
weight=np.array([3,2,1])
ts=np.arange(1979,1990)
coal=WeightedMovingAverage(data,weight,time_series=ts)
coal.show()

class TrendMovingAverage(TimingModel):
    def getTimeSeries(self,ts):
        if type(ts)!=type(None):
            return ts
        else:
            return np.arange(0,self.t+2)
    def __init__(self,arr:np.ndarray,N,time_series=None) -> None:
        super().__init__(arr,time_series)
        self.N=N                                            # 移动平均的项数
        self.time_series=self.getTimeSeries(time_series)    # 时间序列
        self.M1=SimpleMovingAverage(self.arr,self.N).M
        self.M2=SimpleMovingAverage(self.M1,self.N).M
        self.M1t=self.M1[-1:][0]
        self.M2t=self.M2[-1:][0]
        self.a=2*self.M1t-self.M2t
        self.b=2*(self.M1t-self.M2t)/(self.N-1)
    def forecast(self,m=None,target=None):                      
        if type(target)!=type(None):
            this_m=target-int(self.time_series[len(self.time_series)-3])
        elif type(m)!=type(None):
            this_m=m
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return self.a+self.b*this_m
    def show(self,lite=False): 
        print(f'N={self.N}的预测结果如下：')
        raw_data=np.insert(self.arr,self.t,[np.nan,np.nan]).reshape(1,self.t+2)
        est1_data=np.insert(self.M1,self.M1.size,np.nan).reshape(1,self.t+2)
        est2_data=self.M2.reshape(1,self.t+2)
        columns=self.time_series
        index=['观测值','一次预测值','二次预测值']
        data=np.r_[raw_data,est1_data,est2_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'一次移动平均值={round(self.M1t,4)}，二次移动平均值={round(self.M2t,4)}')
        print(f'预测模型为y_(t+m)={round(self.a,4)}+{round(self.b,4)}m')
        print(f'修正预测结果为:')
        print(f'    {self.time_series[len(self.time_series)-2]}={round(self.a+self.b,4)}')
        print(f'    {self.time_series[len(self.time_series)-1]}={round(self.a+2*self.b,4)}')
        print('可以使用forecast()方法进行进一步预测！')

file='E:\DEV\python\math\data\移动平均.xlsx'
data=pd.read_excel(file)
ts=np.arange(1965,1988)
arr=data['发电总量yt'].values

electricity=TrendMovingAverage(arr,6,time_series=ts)
electricity.show(lite=True)
print(electricity.forecast(target=1988))

class ExponentialSmoothing(TimingModel):
    def __init__(self,arr:np.ndarray,alpha,init,time_series=None) -> None:
        super().__init__(arr,time_series)
        self.alpha=alpha
        self.init=init

class PrimaryExponentialSmoothing(ExponentialSmoothing):
    def primaryExponentialSmoothing(self):
        S=np.array([self.init],dtype=np.float64)
        for t in range(0,self.t):
            temp=S.size
            yt=S[-1:][0]
            S=np.insert(S,S.size,self.alpha*self.arr[t]+(1-self.alpha)*yt)
        return S
    def __init__(self,arr:np.ndarray,alpha,init,time_series=None) -> None:
        super().__init__(arr,alpha,init,time_series)
        self.S=self.primaryExponentialSmoothing()
        self.St=self.S[-1:][0]
        self.avg_err=np.mean(self.S[1:]-self.arr)
    def show(self,lite=False):
        print(f'alpha={self.alpha},init={round(self.init,4)}的预测结果如下：')
        raw_data=self.arr.reshape(1,self.t)
        est_data=self.S[1:].reshape(1,self.t)
        columns=self.time_series[:-1]
        index=['观测值','平滑值']
        data=np.r_[raw_data,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'预测结果为：{round(self.St,4)}，平均误差为{round(self.avg_err,4)}')

data=np.array([50,52,47,51,49,48,51,40,48,52,51,59])
init=(data[0]+data[1])/2
alpha=0.2
ts=np.arange(1976,1989)
sale=PrimaryExponentialSmoothing(data,alpha,init,ts)
sale.show()

class QuadraticExponentialSmoothing(ExponentialSmoothing):
    def getTimeSeries(self,ts):
        if type(ts)!=type(None):
            return ts
        else:
            return np.arange(0,self.t+2)
    def __init__(self,arr:np.ndarray,alpha,init,time_series=None) -> None:
        super().__init__(arr,alpha,init,time_series)
        self.time_series=self.getTimeSeries(time_series)    # 时间序列
        self.S1=PrimaryExponentialSmoothing(self.arr,self.alpha,self.init).S[1:]
        self.S2=PrimaryExponentialSmoothing(self.S1,self.alpha,self.init).S[1:]
        self.S1t=self.S1[-1:][0]
        self.S2t=self.S2[-1:][0]
        self.a=2*self.S1t-self.S2t
        self.b=(self.S1t-self.S2t)*self.alpha/(1-self.alpha)
        self.est=(1+1/(1-self.alpha))*self.S1-(1/(1-self.alpha))*self.S2
        self.MSE=np.sum(np.power(self.arr-self.est,2))/self.t           # 均方误差
        self.MAPE=np.sum(np.abs((self.arr-self.est)/self.arr))/self.t   # 平均绝对误差
    def forecast(self,m=None,target=None):                      
        if type(target)!=type(None):
            this_m=target-int(self.time_series[len(self.time_series)-3])
        elif type(m)!=type(None):
            this_m=m
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return self.a+self.b*this_m   
    def show(self,lite=False): 
        print(f'alpha={self.alpha},init={round(self.init,4)}的预测结果如下：')
        raw_data=self.arr.reshape(1,self.t)
        est1_data=self.S1.reshape(1,self.t)
        est2_data=self.S2.reshape(1,self.t)
        est_data=np.insert(self.est,0,np.nan)[:-1].reshape(1,self.t)
        columns=self.time_series[:-2]
        index=['观测值','一次平滑值','二次平滑值','综合估计值']
        data=np.r_[raw_data,est1_data,est2_data,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'一次平滑值={round(self.S1t,4)}，二次平滑值={round(self.S2t,4)}')
        print(f'预测模型为y_(t+m)={round(self.a,4)}+{round(self.b,4)}m,模型均方误差{round(self.MSE,4)},平均绝对误差{round(self.MAPE*100,4)}%')
        print(f'修正预测结果为:')
        print(f'    {self.time_series[len(self.time_series)-2]}={round(self.forecast(1),4)}')
        print(f'    {self.time_series[len(self.time_series)-1]}={round(self.forecast(2),4)}')
        print('可以使用forecast()方法进行进一步预测！')

file='E:\DEV\python\math\data\移动平均.xlsx'
data=pd.read_excel(file)
ts=np.arange(1965,1988)
arr=data['发电总量yt'].values
alpha=0.3
init=arr[0]
electricity=QuadraticExponentialSmoothing(arr,alpha,init,ts)
electricity.show(lite=True)

class TripleExponentialSmoothing(QuadraticExponentialSmoothing):
    def calEst(self):
        A=(3-3*self.alpha+self.alpha**2)/(1-self.alpha)**2
        B=(3-self.alpha)/(1-self.alpha)**2
        C=1/(1-self.alpha)**2
        return A*self.S1-B*self.S2+C*self.S3
    def __init__(self,arr:np.ndarray,alpha,init,time_series=None) -> None:
        super().__init__(arr,alpha,init,time_series)
        self.S3=PrimaryExponentialSmoothing(self.S2,self.alpha,self.init).S[1:]
        self.S3t=self.S3[-1:][0]
        _alpha=self.alpha/(2*(1-self.alpha)**2)
        self.a=3*self.S1t-3*self.S2t+self.S3t
        self.b=_alpha*((6-5*self.alpha)*self.S1t-2*(5-4*self.alpha)*self.S2t+(4-3*self.alpha)*self.S3t)
        self.c=_alpha*self.alpha*(self.S1t-2*self.S2t+self.S3t)
        self.est=self.calEst()
        self.MSE=np.sum(np.power(self.arr-self.est,2))/self.t           # 均方误差
        self.MAPE=np.sum(np.abs((self.arr-self.est)/self.arr))/self.t   # 平均绝对误差
    def forecast(self,m=None,target=None):                      
        if type(target)!=type(None):
            this_m=target-int(self.time_series[len(self.time_series)-3])
        elif type(m)!=type(None):
            this_m=m
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return self.a+self.b*this_m+self.c*this_m**2
    def show(self,lite=False): 
        print(f'alpha={self.alpha},init={round(self.init,4)}的预测结果如下：')
        raw_data=self.arr.reshape(1,self.t)
        est1_data=self.S1.reshape(1,self.t)
        est2_data=self.S2.reshape(1,self.t)
        est3_data=self.S3.reshape(1,self.t)
        est_data=np.insert(self.est,0,np.nan)[:-1].reshape(1,self.t)
        columns=self.time_series[:-2]
        index=['观测值','一次平滑值','二次平滑值','三次平滑值','综合估计值']
        data=np.r_[raw_data,est1_data,est2_data,est3_data,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'一次平滑值={round(self.S1t,4)}，二次平滑值={round(self.S2t,4)}，三次平滑值={round(self.S3t,4)}')
        print(f'预测模型为y_(t+m)={round(self.a,4)}+{round(self.b,4)}m+{round(self.c,4)}m^2,'
              +f'模型均方误差{round(self.MSE,4)},平均绝对误差{round(self.MAPE*100,4)}%')
        print(f'修正预测结果为:')
        print(f'    {self.time_series[len(self.time_series)-2]}={round(self.forecast(1),4)}')
        print(f'    {self.time_series[len(self.time_series)-1]}={round(self.forecast(2),4)}')
        print('可以使用forecast()方法进行进一步预测！')

data=np.array([20.04,20.06,25.72,34.61,51.77,55.92,80.65,131.11,148.58,162.67,232.26])
ts=np.arange(1978,1991)
alpha=0.3
init=np.mean(data[:3])
investment=TripleExponentialSmoothing(data,alpha,init,ts)
investment.show(lite=True)

class FirstDifferenceExponentialSmoothing(ExponentialSmoothing):
    def calEstDif(self):
        est_dif=np.array([self.dif_arr[0]],dtype=np.float64)
        for t in range(1,self.t):
            res=self.alpha*self.dif_arr[t-1]+(1-self.alpha)*est_dif[t-1]
            est_dif=np.insert(est_dif,est_dif.size,res)
        return est_dif[1:]        
    def __init__(self,arr:np.ndarray,alpha,init,time_series=None) -> None:
        super().__init__(arr,alpha,init,time_series)
        self.dif_arr=self.arr[1:]-self.arr[:-1]     # 差分
        self.est_dif=self.calEstDif()               # 差分指数平滑值
        self.est=self.arr[1:]+self.est_dif          # 预测值
        self.est_vel=self.est[-1:][0]
        self.MSE=np.mean(np.power(self.est[:-1]-self.arr[2:],2))                    # 均方误差
        self.MAPE=np.mean(np.abs((self.arr[2:]-self.est[:-1])/self.arr[2:]))        # 平均绝对误差
    def show(self,lite=False):
        print(f'alpha={self.alpha},init={round(self.init,4)}的预测结果如下：')
        raw_data=self.arr.reshape(1,self.t)
        dif_arr=np.insert(self.dif_arr,0,np.nan).reshape(1,self.t)
        est_dif=np.insert(self.est_dif[:-1],0,[np.nan,np.nan]).reshape(1,self.t)
        est_data=np.insert(self.est[:-1],0,[np.nan,np.nan]).reshape(1,self.t)
        columns=self.time_series[:-1]
        index=['观测值','差分','差分平滑值','预测值']
        data=np.r_[raw_data,dif_arr,est_dif,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'预测结果为：{round(self.est_vel,4)}，均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')

data=np.array([24,26,27,30,32,33,36,40,41,44])
ts=np.arange(1977,1988)
alpha=0.4
init=data[0]
fuel=FirstDifferenceExponentialSmoothing(data,alpha,init,ts)
fuel.show()

class AdaptiveFilter(TimingModel):
    def adjWeight(self,w,t,e) -> np.ndarray:
        return w+2*self.k*e*np.flip(self.arr[t-self.N:t])
    def adaptiveFilter(self,_w:np.ndarray,last_sum=0):        
        w=_w
        est_arr=np.array([],dtype=np.float64)   # 预测数列
        e_arr=np.array([],dtype=np.float64)     # 误差数列
        for t in range(self.N,self.t):
            est=np.sum(w*np.flip(self.arr[t-self.N:t]))
            e=self.arr[t]-est
            est_arr=np.insert(est_arr,est_arr.size,est)
            e_arr=np.insert(e_arr,e_arr.size,e)
            w=self.adjWeight(w,t,e)
        # 递归，重复实验减小误差
        if abs(np.sum(e_arr)-last_sum)<1e-9:
            return [est_arr,e_arr,w]
        else:
            return self.adaptiveFilter(w,np.sum(e_arr))
    def __init__(self,arr:np.ndarray,w,k,time_series=None) -> None:
        super().__init__(arr,time_series)
        self.w=np.array(w)                                      # 观测值权重
        self.N=self.w.size                                      # 权数个数
        self.k=k                                                # 学习常数
        self.est_arr=self.adaptiveFilter(self.w)[0]             # 预测数列
        self.e_arr=self.adaptiveFilter(self.w)[1]               # 预测误差
        self.best_w=self.adaptiveFilter(self.w)[2]              # 最佳权重
        self.est_val=np.sum(self.best_w*np.flip(self.arr[-2:])) # 预测值
        self.MSE=np.mean(np.power(self.est_arr-self.arr[self.N:],2))                    # 均方误差
        self.MAPE=np.mean(np.abs((self.arr[self.N:]-self.est_arr)/self.arr[self.N:]))   # 平均绝对误差
    def show(self,lite=False):
        print(f'初始权重w={self.w},学习常数k={self.k}的预测结果如下：')
        raw_data=self.arr.reshape(1,self.t)
        est_data=np.insert(self.est_arr,0,np.full(self.N,np.nan)).reshape(1,self.t)
        columns=self.time_series[:-1]
        index=['观测值','预测值']
        data=np.r_[raw_data,est_data]
        if len(columns)<15 or lite==True:
            display(pd.DataFrame(data,columns=columns,index=index).round(4))
        else:
            display(pd.DataFrame(data.T,columns=index,index=columns).round(4))
        print(f'修正后权重为：{self.best_w.round(4)}')
        print(f'预测结果为：{round(self.est_val,4)}，均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')

data=np.arange(0.1,1.1,0.1)
w=(0.5,0.5)
k=0.9
test=AdaptiveFilter(data,w,k)
test.show()

# 趋势外推
class TrendExtrapolation(TimingModel):
    # 三和法
    def doThreeSum(self,arr):
        S1=np.sum(arr[0:self.m])
        S2=np.sum(arr[self.m:2*self.m])
        S3=np.sum(arr[2*self.m:3*self.m])
        return{'S1':S1,'S2':S2,'S3':S3}
    # 计算系数
    def calCoefficient(self,S1,S2,S3):
        b=np.power((S3-S2)/(S2-S1),1/self.m)
        a=(S2-S1)*(b-1)/(b*(b**self.m-1)**2)
        K=(S1-a*b*(b**self.m-1)/(b-1))/self.m
        return{'b':b,'a':a,'K':K}
    def __init__(self,arr:np.ndarray,time_series=None) -> None:
        super().__init__(arr,time_series)
        if self.arr.size%3!=0:
            raise ValueError('输入观测值数据的长度必须为3的整数倍！')
        self.n=self.arr.size        # 观测值长度
        self.m=int(self.n/3)        # 分期长度
    # 完整预测
    def full_forecast(self):
        res=np.array([],dtype=np.float64)
        for t in range(0,self.t):
            res=np.insert(res,res.size,self.forecast(t+1))
        return res
    # 绘图
    def paint(self):
        plt.figure(figsize=(20,8),dpi=80)
        x=self.time_series[:-1]
        plt.plot(x,self.arr,x,self.est_arr,marker='o')
        plt.legend(labels=["Real Value","Estimated Value"],loc="lower right")
        plt.show()

class ModifiedIndexCurve(TrendExtrapolation):
    def calCascadeGrowthRate(self):
        return (self.arr[2:]-self.arr[1:-1])/(self.arr[1:-1]-self.arr[:-2])
    def forecast(self,t=None,target=None):                      
        if type(target)!=type(None):
            this_t=target-int(self.time_series[0])+1
        elif type(t)!=type(None):
            this_t=t
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return self.K+self.a*self.b**(this_t)
    def __init__(self,arr:np.ndarray,time_series=None) -> None:
        super().__init__(arr,time_series)
        # 逐期增长量比率
        self.cascade_growth_rate=self.calCascadeGrowthRate()
        # 三和法
        self.S1=self.doThreeSum(self.arr)['S1']
        self.S2=self.doThreeSum(self.arr)['S2']
        self.S3=self.doThreeSum(self.arr)['S3']
        # 系数计算
        self.b=self.calCoefficient(self.S1,self.S2,self.S3)['b']
        self.a=self.calCoefficient(self.S1,self.S2,self.S3)['a']
        self.K=self.calCoefficient(self.S1,self.S2,self.S3)['K']
        # 预测
        self.est_arr=self.full_forecast()                                   # 预测数列
        self.MSE=np.mean(np.power(self.est_arr-self.arr,2))                 # 均方误差
        self.MAPE=np.mean(np.abs((self.arr-self.est_arr)/self.arr))         # 平均绝对误差
    def show(self,target=None):
        print('预测结果如下：')
        print('观测值逐期增长量比值落在'+
              f'[{round(np.min(self.cascade_growth_rate),4)},{round(np.max(self.cascade_growth_rate),4)}]范围内;'+
              f'极差为{round(np.ptp(self.cascade_growth_rate),4)},均值为{round(np.mean(self.cascade_growth_rate),4)}'
              )
        print(f'三和法计算结果：S1={round(self.S1,4)},S2={round(self.S2,4)},S3={round(self.S3,4)}')
        print(f'系数计算结果：b={round(self.b,4)},a={round(self.a,4)},K={round(self.K,4)}')
        if type(target)!=type(None):
            print(f'预测结果为：{round(self.forecast(target=target),4)}，均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        else:
            print(f'均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        print(f'可以用forecast()方法进行进一步预测！')
        self.paint()

data=np.array([42.1,47.5,52.7,57.7,62.5,67.1,71.5,75.7,79.8,83.7,87.5,91.1,94.6,97.9,101.1])
ts=np.arange(1969,1985)
radio=ModifiedIndexCurve(data,ts)
radio.show(1986)

class GompertzCurve(TrendExtrapolation):
    def calLogCascadeGrowthRate(self):
        return (np.log(self.arr[2:])-np.log(self.arr[1:-1]))/(np.log(self.arr[1:-1])-np.log(self.arr[:-2]))
    def forecast(self,t=None,target=None):                      
        if type(target)!=type(None):
            this_t=target-int(self.time_series[0])+1
        elif type(t)!=type(None):
            this_t=t
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return self.K*self.a**(self.b**this_t)
    def __init__(self,arr:np.ndarray,time_series=None) -> None:
        super().__init__(arr,time_series)
        # 对数逐期增长量比率
        self.log_cascade_growth_rate=self.calLogCascadeGrowthRate() 
        # 三和法
        self.S1=self.doThreeSum(np.log(self.arr))['S1']
        self.S2=self.doThreeSum(np.log(self.arr))['S2']
        self.S3=self.doThreeSum(np.log(self.arr))['S3']
        # 系数计算
        self.b=self.calCoefficient(self.S1,self.S2,self.S3)['b']
        self.a_=self.calCoefficient(self.S1,self.S2,self.S3)['a']
        self.K_=self.calCoefficient(self.S1,self.S2,self.S3)['K']
        self.a=np.exp(self.a_)
        self.K=np.exp(self.K_)
        # 预测
        self.est_arr=self.full_forecast()                                   # 预测数列
        self.MSE=np.mean(np.power(self.est_arr-self.arr,2))                 # 均方误差
        self.MAPE=np.mean(np.abs((self.arr-self.est_arr)/self.arr))         # 平均绝对误差
    def show(self,target=None):
        print('预测结果如下：')
        print('观测值对数逐期增长量比值落在'+
              f'[{round(np.min(self.log_cascade_growth_rate),4)},{round(np.max(self.log_cascade_growth_rate),4)}]范围内;'+
              f'极差为{round(np.ptp(self.log_cascade_growth_rate),4)},均值为{round(np.mean(self.log_cascade_growth_rate),4)}'
              )
        print(f'三和法计算结果：S1={round(self.S1,4)},S2={round(self.S2,4)},S3={round(self.S3,4)}')
        print(f'系数计算结果：b={round(self.b,4)},a\'={round(self.a_,4)},K\'={round(self.K_,4)},'+
              f'a={round(self.a,4)},K={round(self.K,4)}')
        if type(target)!=type(None):
            print(f'预测结果为：{round(self.forecast(target=target),4)}，均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        else:
            print(f'均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        print(f'可以用forecast()方法进行进一步预测！')
        self.paint()

data=np.array([42.1,47.5,52.7,57.7,62.5,67.1,71.5,75.7,79.8,83.7,87.5,91.1,94.6,97.9,101.1])
ts=np.arange(1969,1985)
radio=GompertzCurve(data,ts)
radio.show(1986)

class LogisticCurve(TrendExtrapolation):
    def calInverseCascadeGrowthRate(self):
        return (1/self.arr[2:]-1/self.arr[1:-1])/(1/self.arr[1:-1]-1/self.arr[:-2])
    def forecast(self,t=None,target=None):                      
        if type(target)!=type(None):
            this_t=target-int(self.time_series[0])+1
        elif type(t)!=type(None):
            this_t=t
        else:
            raise ValueError('target和m至少有一个参数传入！')
        return 1/(self.K+self.a*self.b**this_t)
    def __init__(self,arr:np.ndarray,time_series=None) -> None:
        super().__init__(arr, time_series)
        # 倒数逐期增长量比率
        self.inverse_cascade_growth_rate=self.calInverseCascadeGrowthRate()
        # 三和法
        self.S1=self.doThreeSum(1/self.arr)['S1']
        self.S2=self.doThreeSum(1/self.arr)['S2']
        self.S3=self.doThreeSum(1/self.arr)['S3']
        # 系数计算
        self.b=self.calCoefficient(self.S1,self.S2,self.S3)['b']
        self.a=self.calCoefficient(self.S1,self.S2,self.S3)['a']
        self.K=self.calCoefficient(self.S1,self.S2,self.S3)['K']
        # 预测
        self.est_arr=self.full_forecast()                                   # 预测数列
        self.MSE=np.mean(np.power(self.est_arr-self.arr,2))                 # 均方误差
        self.MAPE=np.mean(np.abs((self.arr-self.est_arr)/self.arr))         # 平均绝对误差
    def show(self,target=None):
        print('预测结果如下：')
        print('观测值倒数逐期增长量比值落在'+
              f'[{round(np.min(self.inverse_cascade_growth_rate),4)},{round(np.max(self.inverse_cascade_growth_rate),4)}]范围内;'+
              f'极差为{round(np.ptp(self.inverse_cascade_growth_rate),4)},均值为{round(np.mean(self.inverse_cascade_growth_rate),4)}'
              )
        print(f'三和法计算结果：S1={round(self.S1,4)},S2={round(self.S2,4)},S3={round(self.S3,4)}')
        print(f'系数计算结果：b={round(self.b,4)},a={round(self.a,4)},K={round(self.K,4)}')
        if type(target)!=type(None):
            print(f'预测结果为：{round(self.forecast(target=target),4)}，均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        else:
            print(f'均方误差为{round(self.MSE,4)}，平均绝对误差为{round(self.MAPE*100,4)}%')
        print(f'可以用forecast()方法进行进一步预测！')
        self.paint()

data=np.array([42.1,47.5,52.7,57.7,62.5,67.1,71.5,75.7,79.8,83.7,87.5,91.1,94.6,97.9,101.1])
ts=np.arange(1969,1985)
radio=LogisticCurve(data,ts)
radio.show(1986)