> CEBG &bjbjَ O"],@RRRRR$6*+"+RRSRRRbR>`*ϻ|Fitting Models to Chaotic Data
J. C. Sprott
Department of Physics, University of Wisconsin, Madison, WI 53706, USA
June 17, 1997
An important problem with many useful applications is finding a mathematical model that replicates certain features of an experimental data record. The model might be useful for prediction, for generation of additional data, for numerical experimentation, or for illuminating the underlying dynamics. The choice of an appropriate model and the method used to fit the model to the data are strongly influenced by the purpose for which the model is to be used. For example, a model optimized for good short-term prediction will usually not do well at capturing the long-term dynamics, and vice versa, especially if the data come from chaotic or noisy dynamics.
This note will describe a numerical algorithm that finds a model equation whose chaotic solutions replicate certain features of the experimental data. The feature chosen is the power spectrum, or more particularly, the auto-correlation function, whose Fourier transform is the power spectrum. Since chaotic data will usually have a broadband power spectrum, the method nearly always produces a model with chaotic solutions, since non-chaotic models would have discrete power spectra.
The model is generally not useful for prediction since it does not preserve the underlying dynamics. In fact, a common technique for generating surrogate data with all the determinism removed is to Fourier transform the data, randomize the phases, and inverse Fourier transform the result, leaving the power spectrum unchanged. However, the method might allow prediction if the model is sufficiently constrained by other factors. It might also be used to exclude certain models by showing that they are incapable of replicating the power spectrum of the data. There may be other applications such as the rapid generation of an unlimited time series with certain spectral properties for testing filters or the bandpass characteristics of communications channels.
The method consists of calculating the auto-correlation function of the data for lags in the range of 1 to 50 (corresponding to frequencies from the Nyquist frequency, 0.5/T, to 0.01/T, where T is the sample time). The same quantity is calculated for the model equation, whose coefficients are then adjusted to give the least square fit. The optimization method is a variant of simulated annealing in which the initial coefficients are randomly chosen, and then neighborhoods in the coefficient space are randomly explored. The size of the neighborhood is initially large but gradually shrinks. Whenever an improved solution is found, its neighborhood is expanded and similarly explored. Both the data and the model are adjusted to have mean 0 and variance 1 before comparisons are made. Initial conditions for the model equation are chosen at 0, and the first 1000 points of the solution are discarded to help ensure that the solution is on the attractor.
Appendix I shows an implementation of the method in BASIC. BASIC was chosen for readability and because when compiled with PowerBASIC for DOS, execution is as fast or faster than with other languages. The program will also run under the much slower Microsoft QuickBASIC. A relatively simple model equation was used of the form
xn+1 = a1 + a2xn + a3xn2 + a4xnxn-1 + a5xn-1 + a6xn-12
This model is easily generalized to higher dimensions (more time lags) or polynomial degree, but it is surprisingly general with appropriate choices of the six coefficients, a1 through a6. The program inputs data from an ASCII file CHAOSFIT.DAT, and it outputs a time series of the same length from the model equations to the file CHAOSFIT.FIT. It will run forever, constantly seeking improved solutions, or until the key is pressed.
An example of the screen output from the program is shown below for a case in which the input data is 2000 points of Gaussian 1/f noise.
The numbers at the top are the optimized values of the six coefficients. The upper trace is the first 50 values of the input data, and the middle trace is the first 50 values of the output data. The two traces at the bottom that almost perfectly overlay are the auto-correlation function of the data and the model solution, respectively. The solutions obtained by this method are usually not unique (other very different values of the coefficients can give similar fits), but they are usually robust to small variations of the coefficients.
Appendix I
'Program CHAOSFIT.BAS fits a chaotic model to experimental data
'(c)1997 by J. C. Sprott (http://sprott.physics.wisc.edu/)
SCREEN 12
jmax% = 50 'Maximum time to correlate
ip% = 50 'Number of time points to plot
kmax% = 6 'Number of coefficients
dmax% = 2 'Number of variables (embedding dimension)
DIM cd(jmax%), cx(jmax%), y(dmax%), a(kmax%), am(kmax%)
OPEN "CHAOSFIT.DAT" FOR INPUT AS #1 'Get length of input data record
imax% = 0
WHILE NOT EOF(1)
imax% = imax% + 1
INPUT #1, d
WEND
CLOSE #1
imax% = imax% - 1
DIM d(imax%), x(imax%)
OPEN "CHAOSFIT.DAT" FOR INPUT AS #1 'Read input data
dav = 0
FOR i% = 0 TO imax%
INPUT #1, d(i%)
dav = dav + d(i%)
NEXT i%
dav = dav / (imax% + 1) 'Average value of input data
CLOSE #1
FOR j% = 0 TO jmax% 'Calculate autocorrelation of input data
cd(j%) = 0
FOR i% = 0 TO imax% - j%
cd(j%) = cd(j%) + (d(i%) - dav) * (d(i% + j%) - dav)
NEXT i%
IF j% = 0 THEN cd0 = cd(0)
cd(j%) = cd(j%) / cd0
NEXT j%
cd0 = SQR(cd0 / (imax% + 1)) 'RMS value of input data
FOR i% = 0 TO imax%: d(i%) = (d(i%) - dav) / cd0: NEXT i%'Normalize to cd0
RANDOMIZE TIMER
rmin = 1E+38: eps = 4: rbest = rmin
WHILE INKEY$ <> CHR$(27) 'Loop until the key is pressed
FOR k% = 1 TO kmax%: a(k%) = am(k%) + eps * (RND - .5): NEXT k%
FOR j% = 1 TO dmax%: y(j%) = 0: NEXT j% 'Zero initial conditions
xav = 0
FOR i% = 1 TO 1000 + imax% 'Generate model equation data
IF ABS(y(1)) > 100 THEN EXIT FOR 'unbounded
y(0) = a(1) + a(2) * y(1) + a(3) * y(1) * y(1) + a(4) * y(1) * y(2) + a(5) * y(2) + a(6) * y(2) * y(2)
IF y(0) = y(1) THEN EXIT FOR 'fixed point
FOR j% = dmax% TO 1 STEP -1: y(j%) = y(j% - 1): NEXT j%
IF i% > 999 THEN x(i% - 1000) = y(1): xav = xav + y(1)
NEXT i%
IF i% = imax% + 1001 THEN 'Calculate autocorrelation of model data
xav = xav / (imax% + 1)
FOR i% = 0 TO imax%: x(i%) = x(i%) - xav: NEXT i% 'Subtract average
FOR j% = 0 TO jmax%
cx(j%) = 0
FOR i% = 0 TO imax% - j%
cx(j%) = cx(j%) + x(i%) * x(i% + j%)
NEXT i%
IF j% = 0 THEN cx0 = cx(0)
cx(j%) = cx(j%) / cx0
NEXT j%
cx0 = SQR(cx0 / (imax% + 1)) 'RMS value of model data
FOR i% = 0 TO imax%: x(i%) = x(i%) / cx0: NEXT i% 'Normalize to cx0
r = 0
FOR j% = 1 TO jmax% 'Calculate sum of squared errors (r)
rsqrt = cd(j%) - cx(j%)
r = r + rsqrt * rsqrt
NEXT j%
IF r <= rmin AND r > 0 THEN 'Test quality of solution
rmin = r
IF rmin < rbest THEN 'A better solution was found
BEEP
rbest = rmin
CLS 'Graph it
LINE (0, 100)-(639, 100), 2
PSET (0, 100 - 40 * d(0))
FOR i% = 1 TO ip%
LINE -(640 * i% / ip%, 100 - 40 * d(i%)), 15
NEXT i%
LINE (0, 240)-(639, 240), 2
PSET (0, 240 - 40 * x(0))
FOR i% = 1 TO ip%
LINE -(640 * i% / ip%, 240 - 40 * x(i%)), 12
NEXT i%
LINE (0, 400)-(639, 400), 2
PSET (0, 300)
FOR j% = 1 TO jmax%
LINE -(640 * j% / jmax%, 400 - 100 * cd(j%)), 15
NEXT j%
PSET (0, 300)
FOR j% = 1 TO jmax%
LINE -(640 * j% / jmax%, 400 - 100 * cx(j%)), 12
NEXT j%
LOCATE 1, 1
FOR k% = 1 TO kmax%: PRINT a(k%), ; : NEXT k%
LINE (0, 0)-(639, 479), , B
END IF
FOR k% = 1 TO kmax%: am(k%) = a(k%): NEXT k% 'Move coefficients
IF eps < 1 THEN eps = 4 * eps
OPEN "CHAOSFIT.FIT" FOR OUTPUT AS #1 'Save model data
FOR i% = 0 TO imax%: PRINT #1, x(i%): NEXT i%
CLOSE #1
ELSE
eps = eps / 1.001 'Shrink the search neighborhood (gradually)
IF eps < .000001 THEN 'Neighborhood too small, start over
rmin = 1E+38
eps = 4
FOR k% = 1 TO kmax%: am(k%) = 0: NEXT k%
END IF
END IF
END IF
WEND
- "#$%&)*+,-.12345679<=>?@BEFGHIKM &&&CJCJ$ jUCJH*CJH*
6CJH*6CJCJ5CJ$B -t MN $$$$ -t MN =Gn,u&.BVlt
'dp5Ei7C'`4T/Cj+9~ + _=Gn,u&.BVlt
''dp5Ei7C'`4T/Cjj+9~ + p !*!L!!!!!"^"v"""")#+ p !*!L!!!!!"^"v"""")#A#_#####:$f$y$$$?%}%%%%3&P&h&&&&&&&&,)#A#_#####:$f$y$$$?%}%%%%3&P&h&&&&&&&&/ =!"#$%nއ9#*XfnPNG
IHDR,gAMA0PLTE{IDATxۖ*E}y$^(Bs;mf@cO'm3=H*"Ce[j[hK eϳvD ތD8Ezq
3ݜ>JrL{{)(A)R$O?2۵HpnWBNvzQ #f2x ~O^&`e?|7n_HPf\z+Z'@J,|~@@V}bcHA@ A
)R HA@ A
)R HA@ A
˶6?`80x
WO4Z )![sa'$^7@Jq 8
Tsa!NM A*_)8 0^_:^iD}/ ConҾ7B
8%_A&u1| JfīcG)`3[t)O-#YOu# pT'`׆D`ݪr|÷;:
J4د\Օ˘͝y}j@)$75u;+^+f@]4"嵈 -**
@&p_J5 o60p)a3T/HSYl-*V/N$:$Ք >;TS4zR
0KiPkge{?ҒNT^))`R J!ޞ9Re@0켦\t!`,GoCMZe)? -6Xh`">p$ZKs# Ff 4Oex۔g<2rgޭ,ǤZ vQMhQvRMM9A32
жR6t&
^GEniP}ڮJ֟"3DJ,-`:m48+j) v
hOvOt0^l%vk=lO9PMh `CpӏE2x}+SZD6m*&ˠG0V7BzZ ^YHgULRxat5<)nRe]8O.ҕKjJh6re½1cא>M۬L{7S %`f/vվr.OӳVg.ݿXZ< 2Sn>=;Z*tq9=pEQ}a]Qu]ܕ$v,vApWX:"&40pj[ Xt?CN XvA@CUTy_^tf\;z**~~CpQN(`ecpqu''N9$<7upN)
*Omu0?,>3G͂F{q@wgu?4{,~vu SzǕ`C `7/=6/6'0ym*8dv_Z`fISrSrΒ.9<31Ge>Cݦٍ8#sP'ϭLTeHҿ3X/ `? 8S<̅Ktxd8EƑ҅kNqD\rWwk2|YLj$:)slr曉'-w`ӐOCǀBGC^Ey|ytN{Kڭ+}og*4Ey1(R|tٰ)XL[[-ɾ8L~N#KPħ:Vc.;XitR/kulLgK9O!dTF-mX hM͂q{@}8V[7FԘ9|04mj$`Ӑ
t
{PDonF
w/IM͂q{@t=h6?g=d/@a\[)H
}8C 3OO0?ʕ^Xpj{
Gpp<o+3sN(8{-SQsK`|2U-ߛ?*[
xA`w'R\`eCxQp=ƱI
X836
ȖcݍV%9V^Tm %? _c_%QFԶWҵ'p섄XR6f\P|1vfF@q5҃3M
_0oOwW0"]F`Z$?u5ii0+90ZAy)jh?>W]~gpD-+Ud6[QaJn*YX2BᶚXSБF^Yo,iK1<{nVA+U+`iDbˋ?roI_ԭjKzSb9f_l,`[+U+씼s8TElɯVZ̢-YRi?4.!
n<6!ԫ7DBr0K*#s
Q;?(2FkxFj
]Qs]'hDטA]Zp
=|;F4<pGCLk
8g}1Pܩ="!IZqRe''`[Ku.+Յ9+T]y;aAU[+(DN8:ۄ.k#oJK
ASX r6xNs=#ZV4[%j\ɷ,-ܨMd
!D|:HuA HA@ A
)R HA@ A
)R HA@ A?AR HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A
)R HA@ A qbIENDB`
[$@$NormalmH 8@8 Heading 1$$@&6CJ4@4 Heading 2$$@&CJ<A@<Default Paragraph Font&>@&Title$CJ0,J@,Subtitle$CJ@C@@Body Text Indent $CJ0Z@"0
Plain TextOJQJ(U@1( Hyperlink>*B*":&'j)#&+ &l,b$އ9#*Xfn:@t(
<
#AB
S ?"8h8(4GQGLnq%*uz&)<AZ]`ctwz}!/18:JM_b!-EISV\adh ;>\aqv%<?BEINjo
dfvxOT"&:>AFHKruCF
} ]!b!!!!!!!C"G"`"c""""
GKnpuy&)PRZ]tw/1EI
;>24<?qs
dgvx34"&$&SU ' ) ` c n!p!!!C"G"`"c""""J. C. Sprott-C:\DOS\TEMP\AutoRecovery save of chaosfit.asdJ. C. Sprott-C:\DOS\TEMP\AutoRecovery save of chaosfit.asdJ. C. SprottC:\chaosfit\chaosfit.docJ. C. SprottC:\chaosfit\chaosfit.docJ. C. SprottC:\chaosfit\chaosfit.docJ. C. SprottC:\WEB\chaos\chaosfit.docJ. C. SprottC:\WEB\chaos\chaosfit.docJ. C. SprottC:\WEB\chaos\chaosfit.docJ. C. SprottC:\WEB\chaos\chaosfit.docJ. C. SprottC:\WINDOWS\TEMP\chaosfit.doc@FF(ÄQQFF"P@GTimes New Roman5Symbol3&Arial?5 Courier New"qhSSBF =#20dD#Fitting Models to Chaotic DataJ. C. SprottJ. C. SprottOh+'0 (
DP\
ht|Fitting Models to Chaotic Dataiitt
J. C. Sprotts t. CNormal.dott
J. C. Sprotts t2 CMicrosoft Word 8.0h@@kH{@$|@$| ՜.+,D՜.+,d hp
UW-Madison Plasma Physics =D#
Fitting Models to Chaotic DataTitle 6>
_PID_GUIDAN{B5A8E608-E6F7-11D0-BCD1-0000C0B4817F}
!"#$%&')*+,-./013456789;<=>?@ADRoot Entry F@,`| Kϻ|F1Table(WordDocumentOSummaryInformation(2DocumentSummaryInformation8:CompObjj
FMicrosoft Word Document
MSWordDocWord.Document.89q~~