Least Squares Method & Error Estimations¶
FIZ353 - Numerical Analysis | 23/10/2020
Emre S. Tasci emre.tasci@hacettepe.edu.tr
Data & Import¶
The free-fall data we will be using is taken from: D. Horvat & R. Jecmenica, "The Free Fall Experiment" Resonance 21 259-275 (2016) [https://doi.org/10.1007/s12045-016-0321-9\].
Here's the contents of our data file:
In [1]:
import pandas as pd
data1 = pd.read_csv("data/03_FreeFallData.csv")
data1.columns
data1
Out[1]:
We don't need the first row and first column, so let's remove them:
In [2]:
data1.drop(0, inplace=True)
data1.drop(11, inplace=True)
data1.drop(['Vert.dist. y/m'],axis=1, inplace=True)
data1
Out[2]:
In [3]:
data1.loc[2,"0.7"]
Out[3]:
In [4]:
type(data1.loc[2,"0.7"])
Out[4]:
In [5]:
data1.dtypes
Out[5]:
In [6]:
data1 = data1.astype('float')
data1.dtypes
Out[6]:
In [7]:
data1.reset_index(inplace=True)
data1
Out[7]:
In [8]:
#data1.reset_index(inplace=True, drop=True)
data1.drop('index',axis=1, inplace=True)
data1
Out[8]:
Plot¶
In [9]:
import seaborn as sns
sns.set_theme() # To make things appear "more beautiful" 8)
In [10]:
# %matplotlib notebook
In [11]:
data2 = data1.copy()
data2
Out[11]:
In [12]:
plt1 = sns.relplot(data=data2,kind="line",marker="o")
k =plt1.set(xticks=data2.index)
In [13]:
data2.mean()
Out[13]:
In [14]:
data_stats = pd.DataFrame(data2.mean())
data_stats.rename(columns={0:'dmean'}, inplace=True )
data_stats['dvar'] = data2.var()
data_stats['dstd'] = data2.std()
data_stats
Out[14]:
In [15]:
data_stats.dstd # Unbiased
Out[15]:
$$\sigma = \sqrt{\frac{\sum_{i}{\left(x_i - \bar{x}\right)^2}}{N}}\;\text{(Population)}$$
or $$\sigma = \sqrt{\frac{\sum_{i}{\left(x_i - \bar{x}\right)^2}}{N-1}}\;\text{(Sample)}$$
In [16]:
import numpy as np
N = data2.shape[0]
for coll in list(data2.columns):
s_dev = 0
s_mean = data2.loc[:,coll].mean()
#print(s_mean)
for x_i in data2.loc[:,coll]:
# print (x_i)
s_dev += (x_i - s_mean)**2
s_dev = np.sqrt(s_dev/(N))
print("{:s}: {:.6f}".format(coll,s_dev))
In [17]:
data2.std(ddof=0) # Biased
Out[17]:
Average over all the sample deviations should be equal to the deviation of the population!
For more information on this Bessel's Correction, check: http://mathcenter.oxford.emory.edu/site/math117/besselCorrection/
Types of Errors¶
True error ($E_t$)¶
$$E_t = \text{true value} - \text{approximation}$$Absolute error ($\left|E_t\right|$)¶
$$\left|E_t\right| = \left|\text{true value} - \text{approximation}\right|$$True fractional relative error¶
$$\frac{\text{true value} - \text{approximation}}{\text{true value}}$$True percent relative error ($\varepsilon_t$)¶
$$\varepsilon_t = \frac{\text{true value} - \text{approximation}}{\text{true value}}\,100\%$$But what if we don't know the true value?..
Approximate percent relative error ($\varepsilon_a$)¶
$$\varepsilon_a = \frac{\text{present approximation} - \text{previous approximation}}{\text{present approximation}}\,100\%$$Computations are repeated until $\left|\varepsilon_a\right| < \left|\varepsilon_s\right|$ ($\varepsilon_s$ is the satisfactory precision criterion).
The result is correct to at least $n$ significant figures given that: $$\varepsilon_s=\left(0.5\times10^{2-n}\right)\,\%$$
Example:¶
$$e^x = 1+x+\frac{x^2}{2}+\frac{x^3}{3!}+\dots+\frac{x^n}{n!}+\dots$$To estimate $e^{0.5}$ so that the absolute value of the approximate error estimate falls below an error criterion conforming to 3 significant figures, how many terms do you need to include?
Solution:
$\varepsilon_s=(0.5\times 10^{2-3})\,\%$
In [18]:
eps_s = 0.5*10**(2-3)
print("{:.2f}%".format(eps_s))
In [19]:
import numpy as np
eps_s = 0.5*10**(2-3)
x = 0.5
eps_a = 1000
e_prev = 0
no_of_terms = 0
print("{:>3} {:^10}\t{:^7}".format("#","e","E_a"))
while(eps_a > eps_s):
e_calculated = e_prev + x**no_of_terms/np.math.factorial(no_of_terms)
eps_a = np.abs(((e_calculated - e_prev)/e_calculated))
print("{:3d}:{:10.6f}\t{:7.4f}".format(no_of_terms+1,e_calculated,eps_a))
e_prev = e_calculated
no_of_terms += 1
In [20]:
import numpy as np
eps_s = 0.5*10**(2-3)
e_sqrt_real = np.sqrt(np.e)
x = 0.5
eps_a = 1000
e_sqrt_prev = 0
no_of_terms = 0
print("{:>3} {:^13}{:^20}{:^12}".format("#","e_calc","E_a","E_t$"))
while(eps_a > eps_s):
e_sqrt_calculated = e_sqrt_prev + x**no_of_terms/np.math.factorial(no_of_terms)
eps_a = np.abs(((e_sqrt_calculated - e_sqrt_prev)/e_sqrt_calculated))*100
eps_t = np.abs((e_sqrt_real - e_sqrt_calculated)/e_sqrt_real)*100
print("{:3d}:{:10.6f}{:16.5f}{:16.5f}".format(no_of_terms+1,e_sqrt_calculated,eps_a,eps_t))
e_sqrt_prev = e_sqrt_calculated
no_of_terms += 1
How good is a mean?¶
Once again, consider our free fall data:
In [21]:
data2
Out[21]:
In [22]:
data_stats
Out[22]:
scipy.optimize.minimize to the rescue!¶
In [23]:
from scipy.optimize import minimize
def fun_err(m,x):
err = 0
for x_i in x:
err += (x_i - m)**2
err = np.sqrt(err/np.prod(x.shape))
return err
In [24]:
fun_err(data2['0.2'].mean(),data2['0.2'])
Out[24]:
In [25]:
fun_err(data2['0.2'].mean()+1,data2['0.2'])
Out[25]:
In [26]:
minimize(fun_err,data2['0.2'].mean(),args=(data2['0.2']))
Out[26]:
In [27]:
data2['0.2'].mean()
Out[27]:
In [28]:
list(data2.columns)
Out[28]:
In [29]:
data_stats.loc['0.2','dmean']
Out[29]:
In [30]:
print("{:^5}: {:^8} ({:^8})".format("col","min","mean"))
for col in list(data2.columns):
res_min = minimize(fun_err,1,args=(data2[col]))
print("{:^5}: {:8.6f} ({:8.6f})".format(col,float(res_min.x),data_stats.loc[col,'dmean']))
In [31]:
print("{:.6f}".format(data_stats.loc[col,'dmean']))
Couldn't the cost function be better?¶
In [32]:
def fun_err2(m,x):
err = 0
for x_i in x:
err += (x_i - m)**2
#err = np.sqrt(err/np.prod(x.shape))
return err
In [33]:
fun_err2(data2['0.2'].mean(),data2['0.2'])
Out[33]:
In [34]:
minimize(fun_err2,data2['0.2'].mean(),args=(data2['0.2']))
Out[34]:
In [35]:
print("{:^5}: {:^8} ({:^8})".format("col","min","mean"))
for col in list(data2.columns):
res_min = minimize(fun_err2,1,args=(data2[col]))
print("{:^5}: {:8.6f} ({:8.6f})".format(col,float(res_min.x),data_stats.loc[col,'dmean']))
In [36]:
def fun_err3(m,x):
err = 0
for x_i in x:
err += np.abs(x_i - m)
#err = np.sqrt(err/np.prod(x.shape))
return err
In [37]:
fun_err3(data2['0.2'].mean(),data2['0.2'])
Out[37]:
In [38]:
minimize(fun_err3,data2['0.2'].mean(),args=(data2['0.2']))
Out[38]:
In [39]:
data_exp = pd.DataFrame(data_stats.dmean)
In [40]:
data_exp
Out[40]:
In [41]:
def freefall_err(g,y_exp,t):
err = 0
y_theo = 0.5*g*t**2
err = (y_theo - y_exp)**2
return np.sum(err)
In [42]:
y_exp = np.array(list(data_exp.index))
print(y_exp)
In [43]:
y_exp.dtype
Out[43]:
In [44]:
y_exp = np.array(list(data_exp.index),dtype=float)
y_exp.dtype
Out[44]:
In [45]:
print(y_exp)
In [46]:
t = np.array(list(data_exp.dmean[:]))
print(t)
We can do that manually!... (???)¶
In [47]:
freefall_err(9,y_exp,t)
Out[47]:
In [48]:
freefall_err(9.1,y_exp,t)
Out[48]:
In [49]:
for g in np.arange(9,10,0.1):
print("{:5.3f}:{:10.6f}".format(g,freefall_err(g,y_exp,t)))
In [50]:
for g in np.arange(9.7,9.9,0.01):
print("{:5.3f}:{:10.6f}".format(g,freefall_err(g,y_exp,t)))
In [51]:
for g in np.arange(9.79,9.81,0.001):
print("{:5.3f}:{:10.8f}".format(g,freefall_err(g,y_exp,t)))
In [52]:
res_min = minimize(freefall_err,x0=1,args=(y_exp,t))
print(res_min)
In [53]:
import matplotlib.pyplot as plt
plt.plot(t,y_exp,"or")
tt = np.linspace(0,0.7,100)
plt.plot(tt,0.5*res_min.x*tt**2,"-b")
plt.show()
Least Squares (scipy.optimize.least_squares) but not least! 8)¶
In [54]:
from scipy.optimize import least_squares
g0 = 1
def fun_yt(g,t):
return 0.5*g*t**2
t = np.array(list(data_exp.dmean[:]))
y_exp = np.array(list(data_exp.index),dtype=float)
def fun(g):
return fun_yt(g,t) - y_exp
res_ls = least_squares(fun,g0)
print(res_ls)
Various Definitions¶
Sum of the squares of the data residuals ($S_t$)¶
$$S_t = \sum_{i}{\left(y_i - \bar{y}\right)^2}$$($\leftrightarrow$ Standard deviation $\sigma=\sqrt{\frac{S_t}{n-1}}$, variance $sigma^2=\frac{\sum_{i}{\left(y_i-\bar{y}\right)^2}}{n-1}=\frac{\sum_{i}{y_i^2}-\left(\sum_{i}{y_i}\right)^2/n}{n-1}$)
Coefficient of variation¶
$$\text{c.v.} = \frac{s_y}{\bar{y}}\times\,100\%$$Sum of the squares of the estimate residuals ($S_r$)¶
$$S_r = \sum_{i}{e_i^2}=\sum_{i}\left(y_i-a_0-a_1 x_i\right)^2$$Standard error of the estimate: ($s_{y/x}$)¶
$$s_{y/x}= \sqrt{\frac{S_r}{n-2}}$$
(a) $S_t$, (b) $S_r$
(Source: Chapra)
Coefficient of Determination ($r^2$)¶
$$r^2 = \frac{S_t-S_r}{S_t}$$Linearization of Non-linear Relationships¶
It doesn't much matter even if the model we're trying to fit is non-linear. We can simply apply a transformation to form it into a linear one. Here are a couple example for handling non-linear functions:
| Model | Nonlinear | Linearized |
|---|---|---|
| exponential | $y=\alpha_1 e^{\beta_1 x}$ | $\ln{y}=\ln{\alpha_1}+\beta_1 x$ |
| power | $y = \alpha_2 x^{\beta_2}$ | $\log{y}=\log{\alpha_2}+\beta_2\log{x}$ |
| saturation-growth-rate | $y=\alpha_3\frac{x}{\beta_3+x}$ | $\frac{1}{y}=\frac{1}{\alpha_3}+\frac{\beta_3}{\alpha_3}\frac{1}{x}$ |
Source: Chapra
Example: Fit the data to the power model (Chapra, 14.6)
Data:
| $i$ | $x_i$ | $y_i$ |
|---|---|---|
| 1 | 10 | 25 |
| 2 | 20 | 70 |
| 3 | 30 | 380 |
| 4 | 40 | 550 |
| 5 | 50 | 610 |
| 6 | 60 | 1220 |
| 7 | 70 | 830 |
| 8 | 80 | 1450 |
Find the optimum $\alpha$ and $\beta$ for the best fit of $y=\alpha x^\beta$ for the given data.
In [55]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
data1 = pd.DataFrame({'i':np.arange(1,9),'x':np.arange(10,90,10),
'y':[25,70,380,550,610,1220,830,1450]})
data1.set_index('i', inplace=True)
In [56]:
data1
Out[56]:
In [57]:
plt.plot(data1.x,data1.y,"o")
plt.show()
Least-squares proper way:¶
We can convert it such that $$\log{y}=\log{\alpha} + \beta\log{x}$$
and as the least square fit for a linear model given as $y' = a_0 + a_1 x'$ is:
$$a_1 = \frac{n\sum{x_i' y_i'} - \sum{x_i'}\sum{y_i'}}{n\sum{x_i'^2}-\left(\sum{x_i'}\right)^2}$$$$a_0 = \bar{y}' - a_1\bar{x}'$$(For derivations, refer to FIZ219 Lecture Notes #5)
and since $x_i' = \log{x_i},\;y_i' = \log{y_i}$:
In [58]:
n = data1.shape[0]
xp = np.log(data1.x)
yp = np.log(data1.y)
a1 = (n*np.sum(xp*yp)-np.sum(xp)*np.sum(yp)) / (n*np.sum(xp**2) - np.sum(xp)**2)
a0 = np.mean(yp) - a1*np.mean(xp)
print("a0: {:7.4f}\na1: {:7.4f}".format(a0,a1))
as $a_0 = \log{\alpha}\rightarrow \alpha = e^{a_0}$ and $a_1 x' = \beta\log{x}\rightarrow \beta = a_1$
In [59]:
alpha = np.exp(a0)
beta = a1
print("alpha: {:7.4f}\nbeta: {:7.4f}".format(alpha,beta))
In [60]:
def fun(alpha, beta, x):
return alpha*x**beta
In [61]:
xx = np.linspace(0,80,100);
yy = fun(alpha,beta,xx)
plt.plot(data1.x,data1.y,"or",xx,yy,"-b")
plt.show()
Minimizing the error function:¶
In [62]:
def fun(alpha, beta, x):
return alpha*x**beta
x = data1.x
y = data1.y
def err(params):
e = y - fun(params[0],params[1],x)
return np.sum(e**2)
from scipy.optimize import minimize
res = minimize(err,[0.274,1.98])
print(res)
alpha2,beta2 = res.x
In [63]:
xx = np.linspace(0,80,100);
yy2 = fun(alpha2,beta2,xx)
plt.plot(data1.x,data1.y,"or",xx,yy2,"-k")
plt.show()
Which one is better?¶
In [64]:
err_ls = err([alpha,beta])
err_min = err([alpha2,beta2])
print("Least-square sum of squares error: {:10.2f}".format(err_ls))
print(" Minimizer sum of squares error: {:10.2f}".format(err_min))
Let's plot the two side by side:
In [65]:
xx = np.linspace(0,80,100);
yy_ls = fun(alpha,beta,xx)
yy_min = fun(alpha2,beta2,xx)
# Blue for least-squares, Black for minimizer
plt.plot(data1.x,data1.y,"or",xx,yy_ls,"-b",xx,yy_min,"-k")
plt.legend(["data","least-squares","minimizer"])
plt.show()
Linear side of things¶
It is always safe to stay on the linear side of the things, whever possible. So, consider the following data:
In [66]:
data_l = np.array(([[ 1. , 4.9142],
[ 2. , 7.1201],
[ 3. , 8.8456],
[ 4. , 10.8113],
[ 5. , 13.2231]]))
data_l
Out[66]:
Let's fit it to a $f(x) = a + bx$ model:
In [67]:
x = data_l[:,0]
y = data_l[:,1]
def fun(x,a,b):
return a+b*x
In [68]:
import scipy.optimize as opt
In [69]:
# Curve fitting
res_cf = opt.curve_fit(fun,x,y)
print(res_cf)
In [70]:
def err(ab):
return np.sum((y - (ab[0]+ab[1]*x))**2)
In [71]:
# Minimizer
res_min = opt.minimize(err,[0.1,0.1])
print(res_min)
In [72]:
# Least-squares
res_ls = opt.least_squares(err,[0.1,0.1],max_nfev=500)
print(res_ls)
In [73]:
xx = np.linspace(0,6,100)
plt.plot(x,y,"or",xx,res_ls.x[0]+res_ls.x[1] * xx,"-b")
plt.show()
Homework¶
- Come up with a (linear) data-set (of 5 to 10 elements) where its mean isn't it's best fit. Why is the one you have found, behaving differently than the g fits we had for each height?
- Understand the derivation of the coefficients $a_0$ and $a_1$ for the linear least-squares fit.
- Calculate the coefficients for the following data for the model $y = a_0 + a_1 x$:
x y 1.0000 4.9142 2.0000 7.1201 3.0000 8.8456 4.0000 10.8113 5.0000 13.2231
using both the least-squares fit and the minimizer methods. Are they different? - Discuss the reason for the difference of the two methods' assignment of the coefficients (in Chapra's 14.6). Do you really think that the minimizer approach is better? Why?
References¶
- D. Horvat & R. Jecmenica, "The Free Fall Experiment" Resonance 21 259-275 (2016) [https://doi.org/10.1007/s12045-016-0321-9]
- This lecture is heavily benefited from Steven Chapra's Applied Numerical Methods with MATLAB: for Engineers & Scientists.