Performing a line fit to observations [ex913.0]ΒΆ

../../_images/sphx_glr_fig-examp-2nd-fit_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 --

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

# Solving Least Squares Problem
W = np.zeros((3,3))
xi = np.zeros((3))
for i, yi in enumerate(y):
    p_i = [1, i, i**2]
    W = W + np.array(np.outer(p_i, p_i))
    xi = xi + np.array(p_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 + x_hat[2]*(k**2)

# -- 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, 6, '$s_i(x) = a_0 + a_1 i + a_2 i = '+str(x_hat[0])+'+'+str(x_hat[1])+ ' i +' + str(x_hat[1])+ 'i^2 $')
axs.legend(loc=1)
axs.set(xlabel='i')
axs.grid(True)

plt.show()

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

Gallery generated by Sphinx-Gallery