denis-bz/1-Holt-doubleexp-smoothing.md

## 1-Holt-doubleexp-smoothing.md

      
    Raw
  

              1-Holt-doubleexp-smoothing.md
            
          
    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

  
## ecg.py
""" 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
""" 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.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
 # 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-ecg.py
""" 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()
	""" 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()