website/archive/binaries/mac/src/disembl/Tisean_3.0.1/source_f/predict.f

   1 c===========================================================================
   2 c
   3 c   This file is part of TISEAN
   4 c
   5 c   Copyright (c) 1998-2007 Rainer Hegger, Holger Kantz, Thomas Schreiber
   6 c
   7 c   TISEAN is free software; you can redistribute it and/or modify
   8 c   it under the terms of the GNU General Public License as published by
   9 c   the Free Software Foundation; either version 2 of the License, or
  10 c   (at your option) any later version.
  11 c
  12 c   TISEAN is distributed in the hope that it will be useful,
  13 c   but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 c   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  15 c   GNU General Public License for more details.
  16 c
  17 c   You should have received a copy of the GNU General Public License
  18 c   along with TISEAN; if not, write to the Free Software
  19 c   Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA  02110-1301  USA
  20 c
  21 c===========================================================================
  22 c   simple nonlinear prediction, fast neighbour search
  23 c   see  H. Kantz, T. Schreiber, Nonlinear Time Series Analysis, Cambridge
  24 c      University Press (1997,2004)
  25 c   author T. Schreiber (1998)
  26 c===========================================================================
  27       parameter(nx=1000000)
  28       dimension x(nx), y(nx)
  29       character*72 file, fout
  30       data eps/0./, frac/0./, ifc/1/
  31       data iverb/1/
  32
  33       call whatido("prediction with locally constant fits",iverb)
  34       id=imust("d")
  35       m=imust("m")
  36       eps=fcan("r",eps)
  37       frac=fcan("v",frac)
  38       ifc=ican("s",ifc)
  39       nmaxx=ican("l",nx)
  40       nexcl=ican("x",0)
  41       jcol=ican("c",0)
  42       isout=igetout(fout,iverb)
  43       if(eps.eq.0.and.frac.eq.0.) call usage()
  44
  45       do 10 ifi=1,nstrings()
  46          call nthstring(ifi,file)
  47          nmax=nmaxx
  48          call readfile(nmax,x,nexcl,jcol,file,iverb)
  49          if(file.eq."-") file="stdin"
  50          if(isout.eq.1) call addsuff(fout,file,"_pred")
  51          call rms(nmax,x,sc,sd)
  52          if(frac.gt.0) eps=sd*frac
  53          iun=istdout()
  54          if(fout.eq." ") iun=istderr()
  55          write(iun,*) "err: ", fcerror(nmax,x,y,m,id,ifc,eps),
  56      .      " "//file(1:index(file," ")-1)
  57  10      call writefile(nmax,y,fout,iverb)
  58       end
  59
  60       subroutine usage()
  61 c usage message
  62
  63       call whatineed(
  64      .   "-d# -m# [-r# | -v#]"//
  65      .   " [-s# -o outfile -l# -x# -c# -V# -h] file(s)")
  66       call ptext("either -r or -v must be present")
  67       call popt("d","delay")
  68       call popt("m","embedding dimension")
  69       call popt("r","absolute radius of neighbourhoods")
  70       call popt("v","same as fraction of standard deviation")
  71       call popt("s","time steps ahead forecast (one step)")
  72       call popt("l","number of values to be read (all)")
  73       call popt("x","number of values to be skipped (0)")
  74       call popt("c","column to be read (1 or file,#)")
  75       call pout("file_pred")
  76       call pall()
  77       stop
  78       end
  79
  80       function fcerror(nmax,y,yp,m,id,ifc,eps)
  81       parameter(im=100,ii=100000000,nx=1000000)
  82       dimension y(nmax),yp(nx),jh(0:im*im),jpntr(nx),nlist(nx)
  83
  84       if(nmax.gt.nx) stop "fcerror: make nx larger."
  85       call base(nmax-ifc,y,id,m,jh,jpntr,eps)
  86       fcerror=0
  87
  88       call rms(nmax,y,sx,sd)
  89       do 10 n=1,(m-1)*id+ifc
  90  10      yp(n)=sx
  91       do 20 n=(m-1)*id+1,nmax-ifc
  92          call neigh(nmax,y,y,n,nmax,id,m,jh,jpntr,eps,nlist,nfound)
  93          av=0
  94          do 30 nn=1,nfound
  95  30         if(nlist(nn).ne.n) av=av+y(nlist(nn)+ifc)
  96          if(nfound.gt.1) then
  97             yp(n+ifc)=av/(nfound-1)
  98          else
  99             yp(n+ifc)=sx
 100          endif
 101  20      fcerror=fcerror+(y(n+ifc)-yp(n+ifc))**2
 102       fcerror=sqrt(fcerror/(nmax-ifc-(m-1)*id))
 103       end
 104