Linear regression - Advanced model¶

Note: this notebook folows the “Linear regression - Linear model” notebook¶

Code author: Emile Roux emile.roux@univ-smb.fr

RISE Slideshow

Scope¶

A finite number \(N\) of data points are available: find the best fit of a given parametrique function going trouth this \(N\) points.

#setup
%load_ext autoreload
%matplotlib nbagg
%autoreload 2
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib as mpl

The data set : a synthetique data set with some noise¶

Nb_data_point = 200 # number of points
xmin, xmax = -.5, 5.5 # x range

x = np.linspace(xmin, xmax, Nb_data_point)
x = x + .2*np.random.rand(Nb_data_point) # add noise

y = x**4-12*x**3+47*x**2-60*x
y = y + 5*np.random.randn(Nb_data_point) #add noise

fig = plt.figure()
plt.plot(x,y,'.',label = "Data")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>

We would like to fit this data with a linear piecewise function¶

To do so we need first a base, let’s use linear piecewise function function use in the FE methode

def hat(x,xc,support=1):
    ### P1 function
    # in:
    # x:
    # xc: center of the peak
    # support:
    y = np.maximum(1-np.abs((x-xc)/support),0.)
    return y

Bases function settings¶

Nf = 15 # number of bases function

xc_list = np.linspace(xmin,xmax,Nf)
xc_list.shape
support = (xmax-xmin)/(Nf-1)

Draw the hat functions¶

N = 1000
x_hat = np.linspace(xmin, xmax, N)

fig = plt.figure()
for i in range(0,len(xc_list)):
    plt.plot(x_hat,hat(x_hat,xc_list[i],support))
plt.grid()

<IPython.core.display.Javascript object>

Fit the data using this base of functions¶

Determination of the coefficent by regression¶

# Construcion of the X matrix
X = np.zeros((Nf,len(x)))
for i in range(0,len(xc_list)):
    X[i,:] = hat(x, xc_list[i], support)
    #X = np.append(X, [hat(x, xc_list[i], support)], axis=0)
X=X.T

# Construcion of the Y matrix
Y = y.T

# Computation of the least square estimator
beta = np.dot(np.linalg.inv(np.dot(X.T,X)),np.dot(X.T,Y))

print("The fitted coeff:")
print(beta)

The fitted coeff:
[ 38.51811083   2.60777129 -16.72832499 -24.50430657 -22.00372903
 -17.16213858  -7.75423456  -4.58088853   1.0449124   -0.6781618
  -1.6077043   -1.08811772  -2.1511271    0.80034793  13.06491387]

Draw the results¶

N = 100
xi = np.linspace(xmin, xmax, N)

yi = np.zeros(xi.shape)
for i in range(0,len(xc_list)):
    yi = yi + beta[i] * hat(xi,xc_list[i],support)

fig = plt.figure()
plt.plot(x,y,'o',label = "Data")
plt.plot(xi,yi,'r',label = "Regression")
plt.legend()
plt.grid()

<IPython.core.display.Javascript object>