Linear Fit [ex9xx.x]

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 *

# -- TEST SIGNAL --

# Constants
K = 50  # number of samples
k = np.array(range(K))

y = np.arange(0,K)*0.1 + .002*np.arange(0,K)**2 + gen_wgn(K, sigma=1, seed=3141)

W = np.zeros((2,2))
xi = np.zeros((2))
for i, yi in enumerate(y):
    W = W + np.array([ [1, i], [i, i**2] ])
    xi = xi + np.array([1, i])*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))

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()

Fitting line model to observations using ALSSMs [ex904.0] Demonstration for a naive implementation of a moving average filter