*
* Example from section 9.5
*
open data motorcycle.dat
data(org=columns,format=free,skip=1) 1 133 index interval accel
*
* The "psi" parameter is actually related to lambda by lambda=exp(2*psi).
*
nonlin psi lambda=exp(2.*psi)
*
dec frml[rect] af
frml af = ||1.0,interval|0.0,1.0||
dec frml[symm] swf
frml swf = ||lambda*interval**3/3.0|lambda*interval**2/2.0,lambda*interval||
*
dlm(a=af,c=||1.0,0.0||,sv=1.0,sw=swf,y=accel,$
     exact,variance=concentrate,type=smooth,method=bfgs) / xstates vstates
*
*  Get the smoothed estimates
*
set smoothed = xstates(t)(1)
set vsmooth  = vstates(t)(1,1)
*
set upper = smoothed+2*sqrt(%variance*vsmooth)
set lower = smoothed-2*sqrt(%variance*vsmooth)
*
spgraph(footer="Fig 9.10. Motorcycle acceleration data analyzed by a cubic spline",vfields=2)
scatter(style=dots,overlay=line,ovcount=3,ovsame) 4
# index accel
# index smoothed
# index lower / 3
# index upper / 3
*
set stdirreg = (accel-smoothed)/sqrt(%variance)
scatter(style=dots)
# index stdirreg
spgraph(done)