""" Hsurvey.py code to make the data for figure ToyH (Fig 5.3, the relative entropy of a Kalman filter) Invoke with: python code/python/Hsurvey.py data Creates files: HtauS and Hsurvey in the data_dir Hsurvey has information for a response surface plot of -cross entropy as a function of sample interval tau_s and log model measurement noise sigma_epsilon HtauS has information for a 2-d figure that has the following plots of cross entropy vs tau_s: One that follows the response surface down its ridge One set at the known simulated measurement noise level A line calculated by Eqn. 5.2 in the book Copyright (c) 2005 Andrew Fraser This file is part of HMM_DS_Code. HMM_DS_Code is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. HMM_DS_Code is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. """ import math, numpy, numpy.linalg, lorenz, EKF, scipy.special, random, sys LAI = numpy.linalg.inv ERF = scipy.special.erf if len(sys.argv) == 2: data_dir = sys.argv[1]+'/' else: data_dir = '' # Parameters s=10.0 r = 28.0 b = 8.0/3 ts = 0.15 Nt = 5000 Ystep = 10**(-4) DevEpsilon = 1e-4 # Square root of the variance of measurement noise SigmaEpsilon = DevEpsilon**2 DevEta = 1e-6 # Square root of the variance of state noise SigmaEta = DevEta**2 # These nine assignments initialize global names results = None Y = None X = None alphaSig = None alphaM = None aM = None F = None AILLY = 0.0 # Average Incremental Log Likelihood of Y def G_func(x): # Function that returns measurement and derivative return numpy.mat(numpy.eye(1))*x[0,0],numpy.mat([[1.0,0.0,0.0]]) def lorstep(x): # Return mu_x(t+1) and derivative F(t) x = x.A.reshape((3,)) # View as array instead of matrix fc,D = lorenz.Ltan_one(x,s,b,r,ts) return numpy.mat(fc.reshape((3,1))),numpy.mat(D) # View as matrices def filter(): global alphaSig,alphaM,aM,aSigma,F,LogStateError, gammaM,gammaS,AILLY alphaSig = Nt*[None] # Updated state covariances alphaM = Nt*[None] # Updated state means aM = Nt*[None] # Forecast state means aSigma = Nt*[None] # Forecast state covariances gammaM = Nt*[None] # Forecast measurement means gammaS = Nt*[None] # Forecast measurement covariances LogStateError = Nt*[None] # Cheat and use the known initial state for the state mean mux = numpy.mat(results[0][0]).T # Initial state variance Sigmax = numpy.mat(numpy.eye(3)*1e-3) # Inverse meas. noise SigEp = numpy.mat(numpy.eye(1))*SigmaEpsilon # Time dependent dynamic noise SigmaEtaT = Nt*[None] for t in xrange(Nt): SigmaEtaT[t] = numpy.eye(3)*SigmaEta # Call the filtering function EKF.ForwardEKF(Y, SigmaEtaT, SigEp, mux, Sigmax, lorstep, G_func, alphaM=alphaM, alphaSig=alphaSig, aM=aM, aSigma=aSigma) # Generate error time series ALE = 0.0 # Average Log X Error AILLY = 0.0 # Average Incremental Log Likelihood of Y Gt = numpy.matrix([1.0,0,0]).T for t in xrange(Nt): xest = aM[t] + alphaM[t] xerr = numpy.mat(X[t]).T - xest fe = numpy.dot(xerr.T,xerr) # Filtered state error squared LogStateError[t] = math.log(fe,10)/2 ALE += LogStateError[t] xt = aM[t][0] sigma_gamma = Gt.T*aSigma[t]*Gt + SigmaEpsilon gammaM[t] = xt gammaS[t] = sigma_gamma #Integrate over bin size here mean = aM[t][0] top = Y[t][0] + Ystep/2 bottom = Y[t][0] - Ystep/2 if sigma_gamma < (Ystep**2)*1e-4: # print 'small sigma_gamma' if bottom < mean and mean < top: like0 = 1.0 else: like0 = 0.0 else: dev = math.sqrt(sigma_gamma*2.0) #2.0 for erf definition z0 = (bottom - mean)/dev z1 = (top - mean)/dev like0 = (ERF(z1) - ERF(z0))/2 z0 = (Y[t][0] - Ystep/2 - mean)/20 z1 = (Y[t][0] + Ystep/2 - mean)/20 like1 = (ERF(z1) - ERF(z0))/2 a = 0.999 # Kludge to limit cost of getting lost like = a*like0 + (1-a)*like1 assert like > 0 AILLY += math.log(like) AILLY/= Nt; def calculate(ev=None): """ Integrate the Lorenz system create X and Y and quantize Y """ global results, X, Y results = [[],[],[],[]] random.seed(3) EKF.TanGen0(results,s=s,r=r,b=b,ts=ts,Nt=Nt, DevEta=DevEta, DevEpsilon=DevEpsilon) X = numpy.array(results[0],numpy.float64) Y = numpy.array(results[1],numpy.float64)/Ystep - 0.5 Y = numpy.ceil(Y)*Ystep # First write intercept and slope for Eqn 5.2 in the book z = 0.5/math.sqrt(2) like = (ERF(z) - ERF(-z))/2 y0 = math.log(like) #The intercept for the "entropy line" is the log #probability of the interval +/- .5 for a normal #distribution. f = open (data_dir+'/HtauS','w') print >>f, "# y = %f -0.906 * t\n" % y0 # Next create data for two lines. One at -4.0 and one on the ridge Tsteps = numpy.arange(.02,.51,.02) dy1 = -4.0 dy2 = -4.85 for ts in Tsteps: # Survey sampling intervals ts DevEpsilon = 1e-10 # Square root of the variance of measurement noise SigmaEpsilon = DevEpsilon**2 calculate() # Create the data DevEpsilon = 10**dy1 SigmaEpsilon = DevEpsilon**2 filter() # Run the EKF to calculate AILLY print >>f, '%4.2f %6.4f'%(ts, AILLY), DevEpsilon = 10**(dy2 + 0.4*ts) SigmaEpsilon = DevEpsilon**2 filter() # Run the EKF to calculate AILLY print >>f, '%6.4f'%AILLY f.close() # Next do a two dimensional survey over sampling interval ts and # DevEpsilon, the scale of the observation noise Tsteps = numpy.arange(.02,.51,.02) dys = numpy.arange(-3.5,-5.6,-.1) f = open(data_dir+'/Hsurvey','w') for ts in Tsteps: DevEpsilon = 1e-10 # Square root of the variance of measurement noise SigmaEpsilon = DevEpsilon**2 calculate() # Make the observations for dy in dys: DevEpsilon = 10**dy SigmaEpsilon = DevEpsilon**2 filter() # Run the EKF to calculate AILLY print >>f, dy, ts, AILLY print >>f # Write an empty line for gnuplot formatting f.close() #Local Variables: #mode:python #End: