Impulse Response Function

In this notebook we further investigate and exemplify the impulse response functions used in time series modelling of groundwater level. The impulse response functions discussed below are used in the GWTSA model.


In [4]:
#Import packages for basic calculations and plotting
import numpy as np
import matplotlib.pyplot as plt
from scipy.special import gammainc, gamma
%matplotlib inline

# First start with the input parameters

t = np.arange(0,3000.,1)
A = 1.
a = 1./10**2.8
n = 1.5

# Create the impulse response curve
X2 = A * ((a ** n) * (t ** (n-1)) * np.exp(-a*t)) / gamma(n)

# Create the step response curve
Fs = A * gammainc(n,t/a)

# Test of response function of Peterson et al. (2014)
A = 1.0
a = 10**2.8
X1 = A*t**(n-1)*np.exp(-t/a)
Fs = np.cumsum(A*t**(n-1)*np.exp(-t/a))

# Create the block response curve
Fb = Fs[1:] - Fs[0:-1]
Fb = np.append(0, Fb) #This is only done for drawing the graph as you normally
# Convolute 

plt.figure(figsize=(10,8))

plt.subplot(321)
plt.bar([0],[1],width=0.1, color='gray')
plt.xlim(-1,8)
plt.ylim(0,1.2)
plt.yticks([1],['$\infty$'])
plt.ylabel('Recharge [L/T]')
plt.subplot(322)
plt.plot(t,X1, 'k')
plt.plot(t,X2, 'r')
plt.ylabel('Response [-]')
plt.xlim(0,3000)


plt.subplot(323)
plt.bar([0],[1],width=8, color='gray')
plt.ylim(0,1.2)
plt.xlim(-1,8)
plt.xticks([0, 8],[0, '$\infty$'])
plt.ylabel('Recharge [L/T]')
plt.subplot(324)
plt.plot(np.arange(len(Fs)),Fs, 'k')
plt.ylabel('Response [-]')
plt.xlim(0,3000)


plt.subplot(325)
plt.bar([0],[1],width=1, color='gray')
plt.ylim(0,1.2)
plt.xlim(-1,8)
plt.ylabel('Recharge [L/T]')
plt.xlabel(' Time [T]')
plt.subplot(326)
plt.plot(np.arange(len(Fb)),Fb, 'k')
plt.ylabel('Response [-]')
plt.xlim(0,3000)
plt.xlabel(' Time [T]')


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-a1efa79ba059> in <module>()
      1 #Import packages for basic calculations and plotting
      2 import numpy as np
----> 3 import matplotlib.pyplot as plt
      4 from scipy.special import gammainc, gamma
      5 get_ipython().magic(u'matplotlib inline')

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in <module>()
   1129 
   1130 # this is the instance used by the matplotlib classes
-> 1131 rcParams = rc_params()
   1132 
   1133 if rcParams['examples.directory']:

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in rc_params(fail_on_error)
    973         return ret
    974 
--> 975     return rc_params_from_file(fname, fail_on_error)
    976 
    977 

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in rc_params_from_file(fname, fail_on_error, use_default_template)
   1098         parameters specified in the file. (Useful for updating dicts.)
   1099     """
-> 1100     config_from_file = _rc_params_in_file(fname, fail_on_error)
   1101 
   1102     if not use_default_template:

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in _rc_params_in_file(fname, fail_on_error)
   1016     cnt = 0
   1017     rc_temp = {}
-> 1018     with _open_file_or_url(fname) as fd:
   1019         try:
   1020             for line in fd:

/Applications/anaconda/lib/python2.7/contextlib.pyc in __enter__(self)
     15     def __enter__(self):
     16         try:
---> 17             return self.gen.next()
     18         except StopIteration:
     19             raise RuntimeError("generator didn't yield")

/Applications/anaconda/lib/python2.7/site-packages/matplotlib/__init__.py in _open_file_or_url(fname)
    998     else:
    999         fname = os.path.expanduser(fname)
-> 1000         encoding = locale.getdefaultlocale()[1]
   1001         if encoding is None:
   1002             encoding = "utf-8"

/Applications/anaconda/lib/python2.7/locale.pyc in getdefaultlocale(envvars)
    541     else:
    542         localename = 'C'
--> 543     return _parse_localename(localename)
    544 
    545 

/Applications/anaconda/lib/python2.7/locale.pyc in _parse_localename(localename)
    473     elif code == 'C':
    474         return None, None
--> 475     raise ValueError, 'unknown locale: %s' % localename
    476 
    477 def _build_localename(localetuple):

ValueError: unknown locale: UTF-8

Interactive example plot

We have now seen a static plot of the impulse, step and block response functions. In the following example an interactive plot is presented. The use can use the sliders to adapt the values of the block response function.


In [3]:
# Define a function that calculates and plots the impulse response function
def impulseresponse_plot(A=1,a=0.01,n=1.5):
    t = np.linspace(0,5000)
    Fs = A * gammainc(n,t/a) #Step Response Function
    Fb = Fs[1:] - Fs[0:-1] # Block Response Function
    Fb = np.append(0,Fb) #This is done for the visualisation, not for actual modelling!!
    plt.plot(t,Fb,'r')
    plt.xlabel('Time [days]')
    plt.ylabel('')
    plt.legend('IRF', loc = 1)
    return Fb

interact(impulseresponse_plot, 
        A = widgets.FloatSliderWidget(min = 0.01, max = 400, step =1, description = 'A', value = 10),
        a = widgets.FloatSliderWidget(min = 1, max = 1000, step =0.005, value = 0.01),
        n = widgets.FloatSliderWidget(min = 0.01, max = 10, step =0.10, value = 1.5))


---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
<ipython-input-3-3a77bc3f755a> in <module>()
     12 
     13 interact(impulseresponse_plot, 
---> 14         A = widgets.FloatSliderWidget(min = 0.01, max = 400, step =1, description = 'A', value = 10),
     15         a = widgets.FloatSliderWidget(min = 1, max = 1000, step =0.005, value = 0.01),
     16         n = widgets.FloatSliderWidget(min = 0.01, max = 10, step =0.10, value = 1.5))

/Applications/anaconda/lib/python2.7/site-packages/IPython/utils/shimmodule.pyc in __getattr__(self, key)
     90             return import_item(name)
     91         except ImportError:
---> 92             raise AttributeError(key)

AttributeError: FloatSliderWidget

In [21]:
t = np.arange(0,3000.,1)
A = 0.3
a = 10**2.8 # 1/a to change for Peterson
n = 1.5
b = 100.

def f(A,a,n,t):
    return A*(t-b)**(n-1)*np.exp(-(t-b)/a) #change /a to * a to comply with Peterson

t_peak = (n-1)*a #time of peak
print t_peak
value_tp = f(A,a,n,t_peak)
print value_tp

def f(A,a,n,t):
    Fs = A * gammainc(n, t/a) # Step response function based on pearsonIII
    Fb = np.append(0,Fs[1:] - Fs[0:-1]) #block reponse function
    #return Fb
    #return (A*t**(n-1)*np.exp(-t/a))/(((n-1)*a)**(n-1)*np.exp(1-n))
    return A*(t-b)**(n-1)*np.exp(-(t-b)/a)

plt.plot(f(A,a,n,t))


315.47867224
3.12973232862
Out[21]:
[<matplotlib.lines.Line2D at 0x10a901790>]

In [5]:
#Create some random shit
from lmfit import minimize, Parameters, Parameter, report_fit
from scipy import random
x = random.rand(len(t))*0.1

f_test = f(A,a,n,t) + x

plt.plot(f_test)


Out[5]:
[<matplotlib.lines.Line2D at 0x1095f5310>]

In [6]:
# Define an objective function
def residuals(X0):
    A = X0['A'].value
    a = X0['a'].value
    n = X0['n'].value
    return (f(A,a,n,t)-f_test)

X0 = Parameters()
#           (Name,  Value,  Vary,   Min,  Max,  Expr)
X0.add_many(('A',   2.0,    True,   None, None,  None),
           ('a',    10**2.8,    True,   None, None,  None),
           ('n',    1.5,    True,   None, None,  None),
)

result = minimize(residuals, X0, method='leastsq')
report_fit(result)

#print residuals(X0)
param = result.params
print param


[[Fit Statistics]]
    # function evals   = 23
    # data points      = 3000
    # variables        = 3
    chi-square         = 2.946
    reduced chi-square = 0.001
[[Variables]]
    A:   1.03716516 +/- 0.001300 (0.13%) (init= 2)
    a:   720.078790 +/- 2.261413 (0.31%) (init= 630.9573)
    n:   1.43336742 +/- 0.002498 (0.17%) (init= 1.5)
[[Correlations]] (unreported correlations are <  0.100)
    C(a, n)                      = -0.879 
    C(A, a)                      = -0.408 
    C(A, n)                      =  0.144 
Parameters([('A', <Parameter 'A', value=1.0371651694278736 +/- 0.0013, bounds=[-inf:inf]>), ('a', <Parameter 'a', value=720.07879040758053 +/- 2.26, bounds=[-inf:inf]>), ('n', <Parameter 'n', value=1.4333674224554911 +/- 0.0025, bounds=[-inf:inf]>)])

In [7]:
x1,x2,x3 = param['A'].value, param['a'].value, param['n'].value
plt.plot(f_test)
plt.plot(f(x1,x2,x3,t), 'r')


Out[7]:
[<matplotlib.lines.Line2D at 0x10959eb50>]

In [8]:
result.covar


Out[8]:
array([[  1.69144846e-06,  -1.20002270e-03,   4.67774435e-07],
       [ -1.20002270e-03,   5.11398882e+00,  -4.96446578e-03],
       [  4.67774435e-07,  -4.96446578e-03,   6.24122484e-06]])

In [ ]:
result.covar

In [ ]: