NN時系列予測

時系列予測: Slide Windows Methos　（C）青山智夫、宮崎大学工学部

slide windowのNNプログラム(digital Fortran v.6)を公開する。

sw.forはこの方法を理解するためのプログラムで、sin(x)の外挿をします。

データも自動生成します。結果はSwtest.xlsです。コメントを読んでください。

Note-1: 階層型NN(multi-layer neural network)による時系列(time-series)データの予測,、slide windows methodについて

この方法の仮定：

1.　時系列データは等時間間隔にサンプリングされている。

2.　過去のデータを集合と考え、その集合を部分集合に分けたとき、その部分集合に未来を予測するデータ成分(関係ではない)が含まれる。

具体的方法：

時系列データを{x0,x1,x2,.....xN}とする。部分集合をたとえば{x0,x1,x2,x3,x4},{x1,x2,x3,x4,x5},{x2,x3,x4,x5,x6},......

とする。次に{x0,x1,x2,x3},{x1,x2,x3,x4},....をNNの入力データ、{x4},{x5},.....を教師データとしてBP学習する。

この学習が可能なのは、部分集合の要素の間の距離が～0でないことである。

学習が完了したNNに{......,xN}を入力して出力される{xN+1}がNNの予測である。

{xN+1}が得られると{....,xN,xN+1}が得られ、それはNNの入力データになるから、さらに{xN+2}が得られる。よってiterativeな予測が可能である。

当然この方法は周期的な現象に適する。実際f(x)=sin(x)やその線形結合関数の予測をするとほとんど完璧な予測ができる。しかも非常に長期の予測が可能である。逆に時間発展現象に対しては威力が発揮できない。

時間発展現象に対しては普通のfunction fittingと同じ近似になる。

増加／減少関係をパターンにして学習すれば(たとえばx0とxiの差分を学習すれば)時間発展現象も予測できるが本来のslide window methodのNNと組合せなければならない。

1 step先だけ予測し次々に測定されている筈の実測データを繰り込んでいく予測を、1 step predictionと言い、NNのslide window method予測にも採用可能である。1 step predictionを短期予測、iterativeにずーと予測する本来の方法を長期予測ともいう。

slide window methodの精度を支配する要因

部分集合{x0,x1,x2,x3,....xM}、window width=M　の決定

データ数N

適切な隠れ層のニューロン数

以上

SW.FORプログラム・ソース　リスト

c This is a program to understand the slide-window method.

c The program includes time-series data ganarator.

C Program: NR4.FOR(version 3), 1995.5.18　このプログラムはPC-9801, MS-FORTRAN 5.1で書かれた

C　それをdigital Fortranにリメイクした。　2000.6.3, 2003.11.29: Tomoo AOYAMA

C NPTN:max.input patterns =< 65

C N60 :max.neurons of the input-layer =< 41

C N25 :max.neurons of the hidden-layer =< 30

C N10 :max.neurons of the output-layer =< 4

C WX :weight matrix between 1st, 2nd layer

C WY :weight matrix between 2nd, 3rd layer

C XIN :rescaled input data Cf. setP routine

C TE :teaching data Cf. gendat, setP routines

C MOPT(1): output list modify cf. PCBK1 routine

c MOPT(20)=NWIN(cf. setP)=slide-window's width

c MOPT(19)=NPTN= number of learning data

c XOPT(3)= termination condition for back-propagation :BPの終了閾値(default:0.0009)

c the main is to set neural network parameters.

IMPLICIT REAL*8 (A-H,O-Z)

COMMON/ARRAYS/XIN(41,65),WX(41,30),H(30),WY(30,4),

1 O(4),DFO(4),DH(30),TE(4,65),TMIN(4),TFAC(4)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

COMMON EX(41)

C base data for random number generation

COMMON/RNDBS/IB,IM,IA,JM,XF

INTEGER*4 IB,IM,IA,JM

IB=173453

IM=37

IA=1

JM=1000000

XF=1./DFLOAT(JM)

OPEN (6,FILE='sw.txt',STATUS='NEW',FORM='FORMATTED')

c default values: eps, epsh are larning coefficients for BP

EPSH=0.2

EPS=EPSH/2.0

DO 77 I=1,20

77 MOPT(I)=0

MOPT(1)=200

c slide-window's width

MOPT(20)=8

c intitial set the number of larning data

c the number is changed as NPTN-"slide window width"

MOPT(19)=21

DO 78 I=1,20

78 XOPT(I)=0.0

XOPT(3)=0.00001

WRITE(6,90)

90 FORMAT(' NR4(v.3)')

c call NN driver

call NNe(XIN,WX,H,WY,O,DFO,DH,TE,TMIN,TFAC)

STOP

END

c NN driver

subroutine NNe(XIN,WX,H,WY,O,DFO,DH,TE,TMIN,TFAC)

implicit real*8 (A-H,O-Z)

dimension XIN(41,65),WX(41,30),H(30),WY(30,4),O(4),

1 DFO(4),DH(30),TE(4,65),TMIN(4),TFAC(4)

c eee= interpolation/extrapolation values

dimension eee(60)

data jj/0/

if(jj.eq.0)then

CALL gendat(TE)

CALL setP(XIN,TE)

CALL LSTP(XIN,TE)

CALL INIT(WX,WY)

c LMTITR= limit number of BP-learning

LMTITR=9001

endif

c BP-learing

CALL PCBK1(LMTITR,XIN,WX,H,WY,O,TE,DFO,DH)

c interpolation (to check original precision of slide-window method)

CALL INTERPOL(eee,XIN,WX,H,WY,O,TE,TMIN,TFAC)

c extrapolation (it is not short-term prediction but long-term)

CALL EXTPOL(eee,XIN,WX,H,WY,O,TE,TMIN,TFAC)

return

end

SUBROUTINE INTERPOL(eee,XIN,WX,H,WY,O,TE,TMIN,TFAC)

c 学習期間再現確認ルーチン

IMPLICIT REAL*8 (A-H,O-Z)

dimension XIN(41,65),WX(41,30),H(30),WY(30,4),O(4),

1 DFO(4),DH(30),TE(4,65),TMIN(4),TFAC(4)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

dimension eee(60)

NWIN=MOPT(20)

N60M=N60-1

WRITE(*,410)

c WRITE(6,410)

410 FORMAT(' interpolation')

DO 100 IPTN=1,NPTN

DO 200 I=1,N25

S=0.

DO 210 J=1,N60

S=S+WX(J,I)*XIN(J,IPTN)

210 CONTINUE

H(I)=1.0/(1.0+DEXP(-S))

200 CONTINUE

DO 220 I=1,N10

S=0.

DO 230 J=1,N25

S=S+WY(J,I)*H(J)

230 CONTINUE

O(I)=S

220 CONTINUE

eee(IPTN)=O(1)

W=eee(IPTN)-TE(2,IPTN+NWIN)

WRITE(*,420)IPTN,eee(IPTN),W

420 FORMAT(1x,I3,':',2G15.6)

WRITE(6,421)IPTN,eee(IPTN),W

421 FORMAT(I3,2G15.6)

c ENDDO

100 CONTINUE

RETURN

END

SUBROUTINE EXTPOL(eee,XIN,WX,H,WY,O,TE,TMIN,TFAC)

c 長期予測ルーチン

IMPLICIT REAL*8 (A-H,O-Z)

dimension eee(60)

dimension XIN(41,65),WX(41,30),H(30),WY(30,4),O(4),

1 DFO(4),DH(30),TE(4,65),TMIN(4),TFAC(4)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

COMMON EX(41)

NWIN=MOPT(20)

N60M=N60-1

WRITE(*,410)

c WRITE(6,410)

410 FORMAT(' extrapolation')

C 外挿のための初期ベクトル作成 to create initial vector for the predictions

DO I=2,N60M

EX(I-1)=XIN(I,NPTN)

ENDDO

EX(N60M)=TE(1,NPTN)

C この要素は定数（バイアス）this is a bias neuron

EX(N60)=1.0

C ここから外挿計算 hereafter, calculations of extrapolation

DO 100 IPTN=1,NWIN*2

DO 200 I=1,N25

S=0.

DO 210 J=1,N60

S=S+WX(J,I)*EX(J)

210 CONTINUE

H(I)=1.0/(1.0+DEXP(-S))

200 CONTINUE

DO 220 I=1,N10

S=0.

DO 230 J=1,N25

S=S+WY(J,I)*H(J)

230 CONTINUE

O(I)=S

220 CONTINUE

eee(IPTN)=O(1)

w=eee(IPTN)-TE(2,NWIN+IPTN +NPTN)

WRITE(*,420)IPTN+NPTN,eee(IPTN),w

420 FORMAT(I3,':',2G15.6)

WRITE(6,421)IPTN+NPTN,eee(IPTN),w

421 FORMAT(I3,2G15.6)

C 外挿用の次のベクトル作成 for next step, creating a vector (: window's data)

DO I=2,N60M

EX(I-1)=EX(I)

ENDDO

EX(N60M)=O(1)

100 CONTINUE

RETURN

END

c wait-routine for debug, nornally it is not used.

subroutine wait0

write(*,90)

90 format(' waitting now')

read(*,100)i

100 format(i1)

return

end

SUBROUTINE gendat(TE)

C data generator: f(x)=(sin(x)+1)/2, 0<x<6pi

c please, to test various time-series phenomena, you change the routine.

IMPLICIT REAL*8 (A-H,O-Z)

dimension TE(4,65)

N=65

PI=0.5*3.141592653589793

DX=PI/(N-1)

DO 100 I=1,N

X=DX*(I-1)

TE(1,I)=(DSIN(X)+1.0)*0.5

100 CONTINUE

RETURN

END

SUBROUTINE LSTP(XIN,TE)

C rescaled input data listing

IMPLICIT REAL*8 (A-H,O-Z)

dimension XIN(41,65),TE(4,65)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

DIMENSION IM(20)

WRITE(*,128)

c WRITE(6,128)

128 FORMAT(' input data')

DO 100 J=1,NPTN

DO 110 I=1,N60

110 IM(I)=9.99*XIN(I,J)

c WRITE(6,130)J,(IM(I),I=1,N60)

WRITE(*,130)J,(IM(I),I=1,N60)

130 FORMAT(1X,I3,2X,49I1)

100 CONTINUE

c WRITE(*,1100)

c READ(*,1110)I

WRITE(*,129)

c WRITE(6,129)

129 FORMAT(' teaching data')

DO 105 J=1,NPTN

c WRITE(6,137)J,(TE(I,J),I=1,N10)

WRITE(*,137)J,(TE(I,J),I=1,N10)

137 FORMAT(1X,I3,2X,10G14.6)

105 CONTINUE

c WRITE(*,1100)

c 1100 FORMAT(' push return key')

c READ(*,1110)I

c 1110 FORMAT(I5)

RETURN

END

SUBROUTINE INIT(WX,WY)

C WX, WY matrix initialize by using random number

IMPLICIT REAL*8 (A-H,O-Z)

dimension WX(41,30),WY(30,4)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

COMMON/RNDBS/IB,IM,IA,JM,XF

INTEGER*4 IB,IM,IA,JM

DO 100 J=1,N25

DO 100 I=1,N60

IB=IB*IM+IA

IB=MOD(IB,JM)

100 WX(I,J)=(XF*DFLOAT(IB)-0.5)*0.2

DO 110 J=1,N10

DO 110 I=1,N25

IB=IB*IM+IA

IB=MOD(IB,JM)

110 WY(I,J)=(XF*DFLOAT(IB)-0.5)*0.2

RETURN

END

SUBROUTINE setP(XIN,TE)

C pattern setting routine, generation of slide-windows

IMPLICIT REAL*8 (A-H,O-Z)

dimension XIN(41,65),TE(4,65)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

NWIN=MOPT(20)

NPTN=MOPT(19)

write(6,306)NPTN,NWIN

306 format(' setP: org.NPTN, NWIN=',2I8)

C window 処理1: ここで入力データを作る creation of "input data" of NN

NPTN=NPTN-NWIN

WRITE(*,301)NPTN

c WRITE(6,301)NPTN

301 FORMAT(' renew.NPTN=',I3)

DO IPTN=1,NPTN

DO I=1,NWIN

XIN(I,IPTN)=TE(1,I+IPTN-1)

ENDDO

C window 処理2: ここで教師データを作る creation of teaching data

C 答え合わせ用データを作る

DO IPTN=1,65

TE(2,IPTN)=TE(1,IPTN)

ENDDO

DO IPTN=1,NPTN

TE(1,IPTN)=TE(2,IPTN+NWIN)

ENDDO

C N10 is the number of output layer neurons. it must be 1

N10=1

C N25 is the number of second layer neurons. it is determined experimentally.

C in case of simple waves, N25 should be 6--8.

c you should check excess-learning on large N25.

N25=6

c N60 is 1st layer's neurons. It includes a bias neuron.

N60=NWIN+1

WRITE(*,430)N60,N25,N10

c WRITE(6,430)N60,N25,N10

430 FORMAT(' N60,N25,N10=',3I8)

DO 410 J=1,NPTN

410 XIN(N60,J)=1.0

900 RETURN

END

SUBROUTINE RSCAL(TE,TMIN,TFAC)

C rescaling of the teaching data

c the routine is not used now.

IMPLICIT REAL*8 (A-H,O-Z)

dimension TE(4,65),TMIN(4),TFAC(4)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

XMIN=10000.

XMAX=-XMIN

I=1

DO 320 J=1,NPTN

IF(TE(I,J).GT.XMAX)XMAX=TE(I,J)

IF(TE(I,J).LT.XMIN)XMIN=TE(I,J)

320 CONTINUE

WRITE(*,321)I,XMAX,XMIN

c WRITE(6,321)I,XMAX,XMIN

321 FORMAT(' RSCAL, DATA',I3,2X,'(MAX,MIN)=',2F12.4)

TFAC(I)=1.0/(XMAX-XMIN)

DO 340 J=1,NPTN

340 TE(I,J)=(TE(I,J)-XMIN)*TFAC(I)

TMIN(I)=XMIN

c WRITE(*,360)TMIN(I)

c WRITE(6,360)TMIN(I)

c 360 FORMAT(' TMIN=',5G13.5)

c WRITE(*,361)TFAC(I)

c WRITE(6,361)TFAC(I)

c 361 FORMAT(' TFAC=',5G13.5)

RETURN

END

SUBROUTINE PCBK1(LMTITR,XIN,WX,H,WY,O,TE,DFO,DH)

C execution of the back propagation

C LMTITR=number of the iteration

c neuron emulate functions are sigmoid (for the 2nd layer) and linear (for output layer)

c it is not normal 3-layer NN.

c This type NN has following characters:

c 1. first learning, 2. extrapolation function, 3. less non-linear-fitting ability, i.e.,

c the NN is hardly trapped by excess-learning

IMPLICIT INTEGER*4 (I-N)

IMPLICIT REAL*8 (A-H,O-Z)

dimension XIN(41,65),WX(41,30),H(30),WY(30,4),O(4),

1 DFO(4),DH(30),TE(4,65)

COMMON/GAB/ETM,XOPT(20)

COMMON/CTL/N60,N25,N10,NPTN,EPS,EPSH,MOPT(20)

DO 90 ITER=1,LMTITR

ERM=-10000.

DO 100 IPTN=1,NPTN

DO 200 I=1,N25

S=0.

DO 210 J=1,N60

S=S+WX(J,I)*XIN(J,IPTN)

210 CONTINUE

H(I)=1.0/(1.0+DEXP(-S))

200 CONTINUE

DO 220 I=1,N10

S=0.

DO 230 J=1,N25

S=S+WY(J,I)*H(J)

230 CONTINUE

C O(I)=1.0/(1.0+DEXP(-S))

O(I)=S

220 CONTINUE

S=0.

DO 320 I=1,N10

W=O(I)-TE(I,IPTN)

C DFO(I)=W*O(I)*(1.0-O(I))

DFO(I)=W

320 S=S+W*W

ERR=S

DO 360 I=1,N25

S=0.

DO 361 J=1,N10

361 S=S+WY(I,J)*DFO(J)

360 DH(I)=S*H(I)*(1.0-H(I))

DO 350 I=1,N10

W=-EPS*DFO(I)

DO 350 J=1,N25

350 WY(J,I)=WY(J,I)+W*H(J)

DO 370 I=1,N25

W=-EPSH*DH(I)

DO 370 J=1,N60

370 WX(J,I)=WX(J,I)+W*XIN(J,IPTN)

IF(ERR.GT.ERM)ERM=ERR

100 CONTINUE

IF(ERM.LT.XOPT(3))go to 9999

90 CONTINUE

WRITE(6,400)ITER,ERM

400 FORMAT(1x,I6,', max.BP.err=',F14.6)

9999 RETURN

END

sin(x)予測結果：

sin(6x), 0<x<2pi, を65点でサンプリングし、1-21点を学習期間、window幅を8とした。

21-8=13が学習期間である。これは予測のための十分なデータが学習期間にある場合である。

1-13までが学習期間、それ以上はNNの予測である。

真値との差は：

学習期間1-13の誤差は+/-0.005以下で、予測になると+/-0.02となり約4倍の誤差となる。

関数fittingによる予測と違い誤差が∞にいかないのが特徴である。

以下、予測のためのデータを減じるとどうなるかを示す。

6pi=3periods, 3*21(学習期間)/65=97%,　予測誤差=0.02, 2%(+16step先)　高精度予測、当然である

5pi=2.5periods, 81%,　0.05, 5%(+16step先)

4pi=2periods, 65%, 0.08 (+6step先)

3pi=1.5, 48%, 0.08 (+6step)

2pi=1.0, 32%, 0.09 (+3step)　NNでも学習期間に将来を予測する成分が無いと、予測不能ということである

0.5pi=0.25, 0.035 (+16step)　関数がf(x)=xに似ているため、予測精度があがる

slide windowは65%とかなり未来を予測するデータが少なくなっても、高精度の予測をします。

その理由は：

SW.FOR　NN予測プログラムは、隠れ層ニューロン関数sigmoid、出力層線形関数です。

NN全体では∑（sigmoid関数）*係数　という近似になってる(*)ので、

原則sigmoid関数によってどのような非線形関数も近似できるはずです。

しかし実際はiteration学習なのでそういかない。

どうしてもsigmoid関数に良く似た関数、sin(x)や線形関数に近似してしまいます。

線形関数に近似するということは、予測を接線方向に延長するということなので、悪くない近似です。

普通人間はそうします。

(*)　このNNは数学的にはRBFと同じです。

学習期間に未来を予測する十分な情報がないとき、現象の性質を使って、情報を補強することもできます。

sin(x)はｙ→-yという変換に対称なので、slide windows+symmetry操作による予測が可能です。

sym操作とは：　{x0,x1,...}と{1-x0,1-x1,....}が対になっていると考え、補データ{1-x0,1-x1,....}を学習すること。

学習期間

2pi=1.0, 32%, 0.09 (+16step)　sym操作により予測精度の向上が見られる

1.5pi=0.75, 24%, 0.09 (+15step)

1pi=0.5, 16%, 0.06 (+16step)

0.5pi=0.25, 0.07 (+16step)　sym操作により予測精度が少し悪くなる

sym操作による予測。F(x)=sin(x), 0<x<2pi(65点),学習期間は1-21点

したがって必要な未来予測情報の32%しかない。

上図の14—がNNによる予測である。誤差は

予測期間14—は次第に増大していくが、十分な精度(10%以下)で予測できている。