Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Holt-Winters double exponential smoothing 2015

A small example of Holt-Winters double exponential smoothing

in Python2 with NumPy and scipy.signal . This is a good bad example, in that yhat[n] == x[n], i.e. predict next == current (a = 1, b = 0) is ~ optimal.

The files:

holt.py:
	ab_BA: a, b -> polynomial coefs B, A
	predict( ab, x ) -> y, predicts x[t+1]
	predict_err( ab, x )

test-holt-ecg.py: runs holt.py on ecg

ecg.py: synthetic ECG data

holt-.9-.01.log

Links:

Holt-Winters smoothing -- nice and clear
Double exponential smoothing

Think about what you really want

absolute vs. relative error
1-step, 2-step ... prediction

Comments welcome, test cases most welcome.

cheers
-- denis

""" http://stackoverflow.com/questions/4387878/simulator-of-realistic-ecg-signal-from-rr-data-for-matlab-or-python """
from __future__ import division
import numpy as np
import scipy.signal as sig
#...............................................................................
rr = np.r_[ 1, 1, .5, 1.5, 1, 1 ] # rr time in seconds
pqrst = sig.wavelets.daub(10) # just to simulate a signal, whatever
# print "wavelets:", pqrst
def ecg( fs=100, rr=rr ):
x = np.concatenate([ sig.resample( pqrst, int(r*fs) ) for r in rr ])
print "ecg: %s * fs %d = %d \n" % (
rr, fs, len(x)) # e.g. [1 1 0.5 1.5 1 1] * fs 100 = 600
return x
""" freqz.py: triv wrap scipy.signal.freqz -> absresponse """
from __future__ import division
import numpy as np
import scipy.signal as sig
# $scipy/signal/signaltools.py lfilter --
#
# a[0]*y[n] = b[0]*x[n] + b[1]*x[n-1] + ... + b[nb]*x[n-nb]
# - a[1]*y[n-1] - ... - a[na]*y[n-na]
#
# (vs R / Autoregressive_model + a[1] ...)
#
# -1 -nb
# b[0] + b[1]z + ... + b[nb] z
# Y(z) = ---------------------------------- X(z)
# -1 -na
# a[0] + a[1]z + ... + a[na] z
#...........................................................................
def freqz( B, A=1, nfreq=9, verbose=1 ):
""" scipy.signal.freqz y = b0 x + b1 x- + b2 x-- ...
- a1 y- ...
-> |response|, print if verbose
"""
B = np.asarray( B )
A = np.asarray( A )
if np.isscalar( nfreq ): # default endpoint=False
T = np.linspace( 0, np.pi, nfreq ) # 6: 0 .1 .. .5
else:
T = nfreq
freqs, response = sig.freqz( B, A, T ) # $scipy/signal/filter_design.py
absresponse = np.abs(response)
if verbose:
print ("freqz: |response| %s \tB %s \tA %s " % (
_ints( absresponse * 100 ), B, A )) .expandtabs( 4 )
# with np.set_printoptions
return absresponse
def _ints( x ):
return np.round( x ).astype(int) # NaN Inf -> - maxint
#...............................................................................
if __name__ == "__main__":
import sys
np.set_printoptions( threshold=100, edgeitems=10, linewidth=100,
formatter = dict( float = lambda x: "%.3g" % x )) # float arrays %.3g
aa = [.1, .2, .5] # iir y += a * (x - y): impulse -> 1 - (1 - a)^n ~ 1/2 at .7/a
nfreq = 9
for arg in sys.argv[1:]:
exec( arg )
for a in aa:
print "\na %.3g --" % a
for b in [ 1,
np.r_[ 1, 1 ] / 2,
np.r_[ 1, 2, 1 ] / 4 ]:
absresponse = freqz( b * a, [1, - (1 - a)], nfreq=nfreq )
print "\nbutter( 2, Wn=.5 ) --"
B, A = sig.butter( 2, Wn=.5 )
freqz( B, A, nfreq=nfreq )
# freqz: |response| [100 100 99 91 71 41 17 4 0] B [0.293 0.586 0.293] A [1 0 0.172]
# # y[n] = .29 x[n] + .59 x[n-1] + .29 x[n-2] - .17 y[n-2]
""" holt.py: Holt-Winters smoothing
ab_BA: a, b -> polynomial coefs B, A
predict( ab, x ) -> y, predicts x[t+1]
predict_err( ab, x )
See
[Holt-Winters smoothing](http://people.duke.edu/~rnau/411avg.htm#HoltLES) -- nice and clear
[Double exponential smoothing](https://en.wikipedia.org/wiki/Exponential_smoothing#Double_exponential_smoothing)
[Hyndman exponentialsmoothing.net](http://exponentialsmoothing.net) -- test csvs
"""
# http://stats.stackexchange.com/questions/6498/seeking-certain-type-of-arima-explanation
from __future__ import division
import sys
import numpy as np
import scipy.signal as sig # lfilter $scipy/signal/signaltools.py
# http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.lfilter.html
__version__ = "2015-04-09 apr denis-bz-py t-online de"
#...............................................................................
def ab_BA( a, b ):
""" -> Holt-Winters a, b -> polynomial coefs B / A (num / denom) for lfilter """
# http://dsp.stackexchange.com/questions/22596/transfer-function-of-double-exponential-smoothing
B = [ a * (1 + b), - a ]
A = [ 1, a * (1 + b) - 2, 1 - a ]
return B, A
def predict( ab, x, zi=None ):
""" -> yhat: Holt-Winters predictor of x[t+1] """
B, A = ab_BA( ab[0], ab[1] )
return sig.lfilter( B, A, x, zi=zi )
def predict_err( ab, x, zi=None, verbose=1 ):
""" rms, yhat, err = predict_err( ab, x ):
yhat = predict( ab, x )
"""
yhat = predict( ab, x, zi=zi )
err = x[1:] - yhat[:-1]
sqrtn = np.sqrt( len(x) - 1 )
rms = np.linalg.norm( err ) / sqrtn
xrms = np.linalg.norm( x[1:] ) / sqrtn
relerr = rms / xrms
if verbose:
print "predict |err| / |x|: %4.1f %% a %-5.3g b %-5.3g xrms %.3g " % (
relerr * 100, ab[0], ab[1], xrms )
return rms, yhat, err
# run-ab: for a ... for b ... test-holt.py
[[ $1 ]] &&
export "$@"
trap exit 2 3 15
for a in .95 .99 # 1 0 best ?
do
for b in .001 .01
do
bold test-holt-ecg.py a=$a b=$b "$@"
log=holt-$a-$b.log
{
From test-holt-ecg.py a=$a b=$b "$@" # pwd, date
#...............................................................................
python test-holt-ecg.py a=$a b=$b "$@"
} > $log 2>&1
egrep -i '^predict|error' $log
done
done
""" test-holt.py """
from __future__ import division
import sys
import numpy as np
import scipy.signal as sig
from freqz import freqz
import holt
import ecg
__version__ = "2015-04-09 apr denis-bz-py t-online de"
#...............................................................................
a, b = .9, .01
fs = 50 # ecg: 6 * fs samples
plot = 0
# to change these params in sh or ipython, run this.py a=1 b=None c=[3] ...
for arg in sys.argv[1:]:
exec( arg )
np.set_printoptions( threshold=100, edgeitems=10, linewidth=120,
formatter = dict( float = lambda x: "%.2g" % x )) # float arrays %.2g
def _ints( x ):
return np.round( x ).astype(int)
#...............................................................................
B, A = holt.ab_BA( a, b )
freqz( B, A ) # B [0.91 -0.9] A [1 -1.1 0.1] -> |response| [100 101 100 ... 83]
x = ecg.ecg( fs ) # <-- your time series here
x *= 100 / x.max()
rms, yhat, err = holt.predict_err( [a, b], x, verbose=1 )
xnext = x[1:]
yhat = yhat[:-1]
title = "Holt-Winters 1-step predictor of synthetic ECG: " \
"|err| %.2g len %s a %.2g b %.3g" % (
rms, len(x), a, b )
print title
print "xnext:", _ints( xnext * 100 )
print "Holt: ", _ints( yhat * 100 )
print "|err|:", _ints( err * 100 )
if plot:
from matplotlib import pyplot as pl, gridspec
import seaborn as sns # nicer plots
sns.set_style("darkgrid")
pl.suptitle( title )
gs = gridspec.GridSpec( 2, 1, height_ratios=[2,1] ) # nrows ncols [left ...]
ax0, ax1 = map( pl.subplot, gs )
# ax0.set_axis_bgcolor( "#f0f0f4" ) # "#EAEAF2" )
# ax1.set_axis_bgcolor( "#f0f0f4" ) # "#EAEAF2" )
t = np.arange(len(x) - 1) / fs
ax0.plot( t, xnext, label="xnext", lw=.5, color="blue" )
ax0.plot( t, yhat, label="Holt", lw=.5, color="orangered" )
ax0.set_xticklabels([])
ax1.set_ylim( -40, 40 )
ax1.vlines( t, 0, err, label="err", color="lightslategrey" )
ax0.legend()
ax1.legend()
if plot >= 2:
from bz.etc import z_numpyutil as nu
nu.savefig( "tmp.png", __file__ )
pl.show()
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.