Fitting line model to observations using ALSSMs [ex904.0]

import numpy as np
import matplotlib.pyplot as plt
from scipy.signal import find_peaks

import lmlib as lm
from lmlib.utils.generator import gen_wgn

# -- TEST SIGNAL --

# Constants
K = 100  # number of samples
k = range(K)

# generate test samples
y = .1*np.arange(0,K) + gen_wgn(K, sigma=.5, seed=3141)


# create ALSSM model with window
alssm_poly = lm.AlssmPoly(poly_degree=1)                     # create ALSSM for polynomial of order 1 (=straight line)
segment = lm.Segment(a=0, b=20, direction=lm.BACKWARD, g=15) # define exponentially decaying window, use backward computations
costseg = lm.CostSegment(alssm_poly, segment)                   # assign window to ALSSM model

se_param = lm.SEParam(costseg)   # prepare memory for intermediate variables and results
se_param.filter(y)  # filter observations y

xs_hat = se_param.minimize_x() # minimize cost function with respect to x



# ----------------  Plot  -----------------
ls_list = ['-','--','-.',':']
ks = [5,45,85] # indeces to be dispayed

#wins = costseg.window(shift=(ks, K))
#trajs = costseg.trajectory(xs_hat[ks], shift=(ks, K), thd=0.1)

trajs = lm.map_trajectories(costseg.trajectories(xs_hat[ks], thd=0.1), ks, K)
wins = lm.map_window(costseg.window(), ks, K)



# axs[0].set( title='noise: $\sigma ='+str(NOISE)+'$')


fig = plt.figure(figsize=(7,4))

axs = fig.add_subplot(5, 1, 1)
for (i, _) in enumerate(wins):
    axs.plot(k, wins[i,:], lw=1, c='k', ls=ls_list[i])
axs.set( ylabel='window(s) \n $w_i$')


axs = fig.add_subplot(5, 1, (2,3))
axs.plot(k, y, 'o', lw=.5, c='k', markersize=2)
for (i, traj) in enumerate(trajs[:,:,0]):
    axs.plot(k, traj, lw=2.0, c='tab:red', ls=ls_list[i], label='$s_i(x_0)$' if i==0 else '')
    axs.scatter(ks[i], traj[ks[i]], marker='x', c='tab:red')
axs.set( ylabel='observations \n $y$')


axs.grid(True)
axs.legend(loc='lower right')
plt.subplots_adjust(hspace=0.4)

axs = fig.add_subplot(5, 1, (4,5))
axs.plot(k, xs_hat[:,0], 'o-', lw=.5, c='tab:green', label='$a_0$ (=offset estimate)', markersize=2)
axs.plot(k, xs_hat[:,1]*100, 'o-', lw=.5, c='tab:blue', label='$a_1 \cdot 100$ (=slope estimate)', markersize=2)
axs.set_ylim((0, 15))
axs.set( ylabel='state estim. \n $x_0 = [a_0, a_1]^T$')
axs.grid(True)
axs.legend(loc='lower right')


plt.show()

Detection of jumps in a signal [ex912.0] Linear Fit [ex9xx.x]