In [1]:
%matplotlib inline
from pylab import *
from numpy import *
from numpy.random import normal, seed
In [2]:
nelem = 40
noise = 0.2

x = linspace(0, 10, nelem)
f = sin(x)
seed(1234)
fn = f + normal(0.0, noise, size=nelem)

subplot(121)
plot(x, fn)

subplot(122)
plot(x[:-1], diff(f))
plot(x[:-1], diff(fn))
Out[2]:
[<matplotlib.lines.Line2D at 0x7fd4140f4358>]
In [3]:
# the antiderivative matrix
def A(n):
    antideriv = ones((n, n))
    antideriv[triu_indices(n, k=1)] = 0.
    return antideriv
In [4]:
print(A(10))
[[ 1.  0.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  0.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  0.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  0.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  0.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  0.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  0.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  0.]
 [ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]]
In [5]:
# check that the inverse of antiderivative approximates derivative
inv(A(10))
Out[5]:
array([[ 1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.]])
In [7]:
AI = A(nelem)
AD = inv(AI)
In [8]:
df = inv(A(nelem)).dot(f)
dfn = inv(A(nelem)).dot(fn)
In [9]:
plot(x, df)
plot(x, dfn)
Out[9]:
[<matplotlib.lines.Line2D at 0x7fd4141b90f0>]
In [10]:
#calculate the condition number
u, s, v = svd(inv(A(nelem)))
k = max(s)/min(s)
print(k)
51.5306511313
In [15]:
# ruzne mu 1. 5.
mu = 5.
xmu = inv(AI.transpose().dot(AI)+mu**2*eye(nelem)).dot(AI.transpose()).dot(fn)
In [16]:
plot(x, xmu)
plot(x, df)
plot(x, dfn)
Out[16]:
[<matplotlib.lines.Line2D at 0x7fd40f471e48>]
In [17]:
mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)
In [18]:
plot(x, xmu)
plot(x, df)
plot(x, dfn)
Out[18]:
[<matplotlib.lines.Line2D at 0x7fd40f3dcbe0>]
In [19]:
from savitzky_golay import *
plot(x, xmu)
plot(x, df)
plot(x, dfn)
plot(x, savitzky_golay(f, 15, 3, 1, rate=1.))
Out[19]:
[<matplotlib.lines.Line2D at 0x7fd40f354828>]
In [ ]:
#another example:
nelem = 40
noise = 0.2

xmin = -1
xmax = 1
x = linspace(xmin, xmax, nelem)
f = abs(x) - abs(xmin)
fn = f + normal(0.0, noise, size=nelem)

AI = A(nelem)
AD = inv(AI)

df = AD.dot(f)
dfn = AD.dot(fn)
In [ ]:
mu = 4.
xmu = inv(AI.transpose().dot(AI)+mu**2*AD.transpose().dot(AD)).dot(AI.transpose()).dot(fn)
In [ ]:
plot(x, savitzky_golay(f, 15, 4, 1, rate=1.))
plot(x, xmu)
plot(x, df)
#plot(x, dfn)
In [ ]: