SE fit of line to samples with ALSSM [ex9xx.x]ΒΆ

Using the line ALSSM model to fit a line model to given samples.

../../_images/sphx_glr_fig-solution-exercise-linear-fit-alssm-simple_001.png
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_slopes, gen_wgn
from lmlib.utils.generator import *

# -- Generating test signal --
K = 50  # number of samples
k = np.arange(0, K)
y = k*0.1 + .002*k**2 + gen_wgn(K, sigma=1, seed=3141)


# --- Definition of an ALSSM for a line model ---
A = np.array([[1, 1], [0, 1]])
c = np.array([[1, 0]]) # row vector

# --- Initializing sum terms ---
N = A.shape[0] # get ALSSM system order (=matrix dimension)
W = np.zeros((N, N))
xi = np.zeros((N,1))
Ai = np.eye(N)  # initialize A^0 as the unity matrix

# --- Perform recursive summing ---
for i, yi in enumerate(y):
    Ai = Ai@A # recursive computation of A^i

    W = W + Ai.transpose()@c.transpose()@c@Ai
    xi = xi + Ai.transpose()@c.transpose()*yi



x_hat = np.linalg.inv(W)@xi # estimate of slope and offset
y_hat =  x_hat[0]*np.ones((K)) + x_hat[1]*k


# ---- PLOTTING ----
_, axs = plt.subplots(1, 1, sharex='all', figsize=(6, 3))
k = np.array(range(K))
axs.plot(k, y_hat, c='r', lw=1.0, label='$\hat{y}_i$')
axs.plot(k, y, 'o', c='b', lw=1.0, label='$y_i$')
axs.text(0, 9, '$s_i(x) = a_0 + a_1 i = '+str("%.4f" % x_hat[0])+'+'+str("%.4f" % x_hat[1])+ 'i $')
axs.legend(loc=4)
axs.set(xlabel='i')
axs.grid(True)
plt.show()

Total running time of the script: ( 0 minutes 0.159 seconds)

Gallery generated by Sphinx-Gallery