前言
从60年代初期以来自回归模型就用于水资源规划和设计中，这类模型之所以在水文学中得到重视并有较大的吸引力，其主要原因是它们具有时间相依的非常直观的形式，同时建立模型和具体应用简单。
模型简介
AR模型是一种线性预测，即已知N个数据，可由模型推出第N点前面或后面的数据（设推出P点），所以其本质类似于插值，其目的都是为了增加有效数据，只是AR模型是由N点递推，而插值是由两点（或少数几点）去推导多点，所以AR模型要比插值方法效果更好。
一般自回归模型通常表示为下式：
模型实现
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 | #!/usr/bin/env python # -*- coding: utf-8 -*- # @Date : 2017-05-05 21:31:16 # @Author : Edward Linus ([email protected]) # @Link : https://edlinus.cn # @Version : 0.1 import os import numpy as np ''' 计算Cvx k:模比系数矩阵 ''' def cvx(k): cv=0 n=len(k);sum=0 tmp=k-1 tmp=np.power(tmp,2) for num in tmp: sum+=num cv=np.sqrt(sum/(n-1)) return cv ''' 计算Csx k:模比系数矩阵 cv:离势系数 ''' def csx(k,cv): cs=0 n=len(k);sum=0 tmp=k-1; tmp=np.power(tmp,3) for num in tmp: sum+=num cs=sum/(np.power(cv,3)*(n-3)) return cs ''' 计算相关系数 x:系列 ss:下标号 ''' def rho(x,ss): xt=x[ss:] xl=x[:len(x)-ss] cor=np.corrcoef(xt,xl) return cor ''' 计算h回归系数 x:系列 n:阶数 ''' def fhi(x,n): r=np.zeros([n]) f=np.zeros([n,n]) for i in range(1,n+1): r[i-1]=rho(x,i)[0,1] f[0,0]=r[0]#初值 for j in range(1,n): cr=cross(f,r,j) f[j,j]=(r[j]-cr[0])/(1-cr[1]) for k in range(0,j): f[j,k]=f[j-1,k]-f[j,j]*f[j-1,j-k-1] return f,r ''' 交叉乘 m1:系列1 m2:系列2 n:位置 ''' def cross(m1,m2,n): sum_1=0;sum_2=0 for i in range(0,n): sum_1+=m1[n-1,i]*m2[n-i-1] sum_2+=m1[n-1,i]*m2[i] return sum_1,sum_2 ''' 计算回归系数和相关系数的积 f:回归系数矩阵 r:相关系数矩阵 n:矩阵长度 ''' def fhirho(f,r,n): sum=0 for i in range(0,n): sum+=f[n-1,i]*r[i] return sum ''' 计算最小值及下标 m:系列 ''' def min(m): i=0;mnum=m[i] for x in range(1,len(m)): if m[x]<mnum: i=x mnum=m[x] return mnum,i ''' 计算数值在数组中的位置 num:数字 list:数组 ''' def pos(num,list): ub = 0;db=0 for i in range(0,len(list)-1): ubn=list[i];dbn=list[i+1] if num>=ubn and num<=dbn : ub=i;db=i+1 return ub,db ''' 生成PIII型分布的随机误差 u:0-1均匀分布的随机数 csr:随机误差的离势系数 plist:φ值表矩阵 ''' def p3num(u,csr,dr,plist): u=u*100 xp=pos(csr,plist[:,0]) yp=pos(u,plist[0,:]) if xp[0]==xp[1] or yp[0]==yp[1]: return 0 cst1=plist[xp[0],yp[0]]+plist[xp[1],yp[0]]*(u-plist[0,yp[0]])/(plist[0,yp[1]]-plist[0,yp[0]]) cst2=plist[xp[0],yp[1]]+plist[xp[1],yp[1]]*(u-plist[0,yp[1]])/(plist[0,yp[1]]-plist[0,yp[0]]) p3n=cst1+(cst2-cst1)*(csr-plist[xp[0],0])/(plist[xp[1],0]-plist[xp[0],0]) p3n=dr*p3n return p3n filepath=os.path.dirname(os.path.realpath(__file__)) x=np.loadtxt(filepath+"\\x.dat")#数据导入 n=len(x) #数据数目 pmax=int(n/5) #最大阶数上界 pmin=int(n/10) #最大阶数下界 fpe=np.zeros([pmax-pmin+1]) #模型误差矩阵 xm=np.mean(x) #平均值 y=x-xm #中心化系列 k=x/xm #模比系数 cv=cvx(k) #离势系数 cs=csx(k,cv) #偏态系数 dx=np.std(x) #标准差 ''' flag=0 for p in range(pmin,pmax+1): (f,r)=fhi(x,p) #回归系数及相关系数矩阵 dk=np.power(dx,2)*(1-fhirho(f,r,len(r))) #均方差 fperr=dk*(1+3*p/n) #λ=3时的模型误差 fpe[flag]=fperr flag+=1 p=min(fpe)[1]+pmin ''' p=pmax #TODO (f,r)=fhi(x,p) cof=f[len(r)-1,:] #回归系数 b=xm #回归方程截距 str1="" for c in range(0,len(cof)): b+=-cof[c]*xm if cof[c]>=0: str1+='+'+str(cof[c]) + "x|t-" + str(c+1) else: str1+=str(cof[c]) + "x|t-" + str(c+1) str0="回归方程为：y=" str0+=str(b) + str1 + "+rt" print(str0) plist=np.loadtxt(filepath+"\\φ.dat") r_2=0 #R^2 for i in range(0,p): r_2+=f[p-1,i]*r[i] dr=dx*np.sqrt(1-r_2) #随机变量的标准差 csr=cs*(1-np.power(r_2,1.5))/np.power(1-r_2,1.5) #随机变量的偏态系数 fp = input("请输入预报长度:\n") #预报长度 #fp=10 yf=np.zeros([fp]) xin=x[len(x)-len(cof):] for i in range(0,fp): u=np.random.rand() p3n=p3num(u,csr,dr,plist) #print(p3n) yf[i]+=b+p3n for j in range(0,len(cof)): yf[i]+=xin[j]*cof[len(cof)-j-1] for k in range(0,len(cof)-1): xin[k]=xin[k+1] xin[len(cof)-1]=yf[i] np.savetxt(filepath+"\\y.dat",yf,fmt='%s') print("预报结果已保存在"+filepath+"\\y.dat") |
插值所需要的φ值数据如下
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | -0.0 0.01 0.10 0.20 0.33 0.50 1.0 2.0 5.0 10.0 20.0 50.0 75.0 90.0 95.0 99.0 0.0 3.72 3.09 2.88 2.71 2.58 2.33 2.05 1.64 1.28 0.84 0.00 -0.67 -1.28 -1.64 -2.33 0.1 3.94 3.23 3.00 2.82 2.67 2.40 2.11 1.67 1.29 0.84 -0.02 -0.68 -1.27 -1.62 -2.25 0.2 4.16 3.38 3.12 2.92 2.76 2.47 2.16 1.70 1.30 0.83 -0.03 -0.69 -1.26 -1.59 -2.18 0.3 4.38 3.52 3.24 3.03 2.86 2.54 2.21 1.73 1.31 0.82 -0.05 -0.70 -1.24 -1.55 -2.10 0.4 4.61 3.67 3.36 3.14 2.95 2.62 2.26 1.75 1.32 0.82 -0.07 -0.71 -1.23 -1.52 -2.03 0.5 4.83 3.81 3.48 3.25 3.04 2.68 2.31 1.77 1.32 0.81 -0.08 -0.71 -1.22 -1.49 -1.96 0.6 5.05 3.96 3.60 3.35 3.13 2.75 2.35 1.80 1.33 0.80 -0.10 -0.72 -1.20 -1.45 -1.88 0.7 5.28 4.10 3.72 3.45 3.22 2.82 2.40 1.82 1.33 0.79 -0.12 -0.72 -1.18 -1.42 -1.81 0.8 5.50 4.24 3.85 3.55 3.31 2.89 2.45 1.84 1.34 0.78 -0.13 -0.73 -1.17 -1.38 -1.74 0.9 5.73 4.39 3.97 3.65 3.40 2.96 2.50 1.86 1.34 0.77 -0.15 -0.73 -1.15 -1.35 -1.66 1.0 5.96 4.53 4.09 3.76 3.49 3.02 2.54 1.88 1.34 0.76 -0.16 -0.73 -1.13 -1.32 -1.59 1.1 6.18 4.67 4.20 3.86 3.58 3.09 2.58 1.89 1.34 0.74 -0.18 -0.74 -1.10 -1.28 -1.52 1.2 6.41 4.81 4.32 3.95 3.66 3.15 2.62 1.91 1.34 0.73 -0.19 -0.74 -1.08 -1.24 -1.45 1.3 6.64 4.95 4.44 4.05 3.74 3.21 2.67 1.92 1.34 0.72 -0.21 -0.74 -1.06 -1.20 -1.38 1.4 6.87 5.09 4.56 4.15 3.83 3.27 2.71 1.94 1.33 0.71 -0.22 -0.73 -1.04 -1.17 -1.32 1.5 7.09 5.23 4.68 4.24 3.91 3.33 2.74 1.95 1.33 0.69 -0.24 -0.73 -1.02 -1.13 -1.26 1.6 7.31 5.37 4.80 4.34 3.99 3.39 2.78 1.96 1.33 0.68 -0.25 -0.73 -0.99 -1.10 -1.20 1.7 7.54 5.50 4.91 4.43 4.07 3.44 2.82 1.97 1.32 0.66 -0.27 -0.72 -0.97 -1.06 -1.14 1.8 7.76 5.64 5.01 4.52 4.15 3.50 28.50 1.98 1.32 0.64 -0.28 -0.72 -0.94 -1.02 -1.09 1.9 7.98 5.77 5.12 4.61 4.23 3.55 2.88 1.99 1.31 0.63 -0.29 -0.72 -0.92 -0.98 -1.04 2.0 8.21 5.91 5.22 4.70 4.30 3.61 2.91 2.00 1.30 0.61 -0.310 -0.710 -0.895 -0.949 -0.989 2.1 8.43 6.04 5.33 4.79 4.37 3.66 2.93 2.00 1.29 0.59 -0.320 -0.710 -0.869 -0.914 -0.945 2.2 8.65 6.17 5.43 4.88 4.44 3.71 2.96 2.00 1.28 0.57 -0.330 -0.700 -0.844 -0.879 -0.905 2.3 8.87 6.30 5.53 4.97 4.51 3.76 2.99 2.00 1.27 0.55 -0.340 -0.690 -0.820 -0.849 -0.867 2.4 9.08 6.42 5.63 5.05 4.58 3.81 3.02 2.01 1.26 0.54 -0.350 -0.680 -0.795 -0.820 -0.831 2.5 9.30 6.55 5.73 5.13 4.65 3.85 3.04 2.01 1.25 0.52 -0.360 -0.670 -0.772 -0.791 -0.800 2.6 9.51 6.67 5.82 5.20 4.72 3.89 3.06 2.01 1.23 0.50 -0.370 -0.660 -0.748 -0.764 -0.769 2.7 9.72 6.79 5.92 5.28 4.78 3.93 3.09 2.01 1.22 0.48 -0.370 -0.650 -0.726 -0.736 -0.740 2.8 9.93 6.91 6.01 5.36 4.84 3.97 3.11 2.01 1.21 0.46 -0.380 -0.640 -0.702 -0.710 -0.714 2.9 10.14 7.03 6.10 5.44 4.90 4.01 3.13 2.01 1.20 0.44 -0.390 -0.630 -0.680 -0.687 -0.690 |
