Hysteresis and regression analysis example¶
Given a calibration of an instrument for an increasing and decreasing input \(q_i\) [kPa] and output of the instrument \(q_o\) [kPa]
import pylab
from scipy.stats import t
import numpy as np
from pylab import array, figure, grid, legend, plot, r_, size, title, xlabel, ylabel
pylab.rcParams['figure.figsize'] = 10, 8
pylab.rcParams['font.size'] = 16
from IPython.core.display import Image
Image(filename='images/pressure_calibration_example.png',width=400)
# increasing
xi = array([0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0])
yi = array([-1.12, 0.21, 1.18, 2.09, 3.33, 4.50, 5.26, 6.59, 7.73, 8.68, 9.80])
# decreasing
xd = xi.copy()
yd = array([-0.69, 0.42, 1.65, 2.48, 3.62, 4.71, 5.87, 6.89, 7.92, 9.10, 10.20])
figure()
plot(xi,yi,'r:o',xd,yd,'b:s')
grid('on')
xlabel('True $P$ [kPa]')
ylabel('Measured $P$ [kPa]')
legend(('increasing','decreasing'))
Estimate the hysteresis error:¶
\(e_h = y_{up} - y_{down}\)
\(e_{h_{max}} = max(|e_h|)\)
$e_{h_{max}}\% = 100\% \cdot \frac{e_{h_{max}}}{y_{max}-y_{min}} $
e_h = [-0.43 -0.21 -0.47 -0.39 -0.29 -0.21 -0.61 -0.3 -0.19 -0.42 -0.4 ] [kPa]
e_hmax= 0.61 [kPa]
e_hmax_p = 100*e_hmax/(np.max(xi) - np.min(xd))
print("Relative error = %3.2f%s FSO" % (e_hmax_p,"%"))
Relative error = 6.10% FSO
Regression analysis¶
Following the recipe of http://www.answermysearches.com/how-to-do-a-simple-linear-regression-in-python/124/
def linreg(X, Y):
"""
Summary
Linear regression of y = ax + b
Usage
real, real, real = linreg(list, list)
Returns coefficients to the regression line "y=ax+b" from x[] and y[], and R^2 Value
"""
N = len(X)
if N != len(Y): raise(ValueError, 'unequal length')
Sx = Sy = Sxx = Syy = Sxy = 0.0
for x, y in zip(X, Y):
Sx = Sx + x
Sy = Sy + y
Sxx = Sxx + x*x
Syy = Syy + y*y
Sxy = Sxy + x*y
det = Sx * Sx - Sxx * N # see the lecture
a,b = (Sy * Sx - Sxy * N)/det, (Sx * Sxy - Sxx * Sy)/det
meanerror = residual = residualx = 0.0
for x, y in zip(X, Y):
meanerror = meanerror + (y - Sy/N)**2
residual = residual + (y - a * x - b)**2
residualx = residualx + (x - Sx/N)**2
RR = 1 - residual/meanerror
# linear regression, a_0, a_1 => m = 1
m = 1
nu = N - (m+1)
sxy = np.sqrt(residual / nu)
# Var_a, Var_b = ss * N / det, ss * Sxx / det
Sa = sxy * np.sqrt(1/residualx)
Sb = sxy * np.sqrt(Sxx/(N*residualx))
# We work with t-distribution, ()
# t_{nu;\alpha/2} = t_{3,95} = 3.18
tvalue = t.ppf(1-(1-0.95)/2, nu)
print("Estimate: y = ax + b")
print("N = %d" % N)
print("Degrees of freedom $\\nu$ = %d " % nu)
print("a = %.2f $\\pm$ %.3f" % (a, tvalue*Sa/np.sqrt(N)))
print("b = %.2f $\\pm$ %.3f" % (b, tvalue*Sb/np.sqrt(N)))
print("R^2 = %.3f" % RR)
print("Syx = %.3f" % sxy)
print("y = %.2f x + %.2f $\\pm$ %.3f V" % (a, b, tvalue*sxy/np.sqrt(N)))
return a, b, RR, sxy
print("Increasing")
ai,bi,Ri,sxyi = linreg(xi,yi)
print(ai,bi,Ri)
print("Decreasing")
ad,bd,Rd,sxd = linreg(xd,yd)
print(ad,bd,Rd)
Increasing Estimate: y = ax + b N = 11 Degrees of freedom $\nu$ = 9 a = 1.08 $\pm$ 0.008 b = -1.03 $\pm$ 0.045 R^2 = 0.999 Syx = 0.116 y = 1.08 x + -1.03 $\pm$ 0.079 V 1.0823636363636362 -1.025454545454545 0.9990647108930198 Decreasing Estimate: y = ax + b N = 11 Degrees of freedom $\nu$ = 9 a = 1.08 $\pm$ 0.005 b = -0.67 $\pm$ 0.028 R^2 = 1.000 Syx = 0.072 y = 1.08 x + -0.67 $\pm$ 0.049 V 1.082272727272727 -0.668636363636363 0.9996380668857191
figure()
plot(xi,yi,'ro',xi,xi*ai+bi,'r--',xd,yd,'bs',xd,xd*ad+bd,'b--')
xlabel('$P_{true}$ [kPa]')
ylabel('$P_{meas}$ [kPa]')
title('Regression')
leg1 = '%3.2f $x$ %3.2f' % (ai,bi)
leg2 = '%3.2f $x$ %3.2f' % (ad,bd)
legend(('incr.',leg1,'decr.',leg2),loc='best')
# all all together:
a,b,R,_ = linreg(r_[xi,xd],r_[yi,yd])
plot(xi,xi*a+b,'k-',lw=2)
leg3 = '$q_o$ = %3.2f $q_i$ %3.2f' % (a,b)
legend(('incr.',leg1,'decr.',leg2,leg3),loc='best')
Estimate: y = ax + b N = 22 Degrees of freedom $\nu$ = 20 a = 1.08 $\pm$ 0.006 b = -0.85 $\pm$ 0.037 R^2 = 0.997 Syx = 0.208 y = 1.08 x + -0.85 $\pm$ 0.093 V
\[q_0 = m q_i + b\]
\[q_0 = 1.8 q_i - 0.85\]
\[\sigma_m^2 = \frac{N \sigma_{q_0}^2}{N\sum q_i^2 - \left(\sum q_i \right)^2}\]
\[\sigma_b^2 = \frac{\sigma_{q_0}^2\sum q_i^2}{N\sum q_i^2 - \left(\sum q_i \right)^2}\]
\[\sigma_{q_0}^2 = \frac{1}{N}\sum\left(m\,q_i + b - q_0 \right)^2 \]
\[\sigma_{q_i}^2 = \frac{1}{N}\sum\left( \frac{q_0 - b}{m} - q_i \right)^2 = \frac{\sigma_{q_0}^2}{m^2} \]
def regression_analysis(qi,qo):
m,b,R,syx = linreg(qi,qo)
# print 'qo = %3.2f qi %3.2f' % (m,b)
N = size(qi)
d = m*qi + b - qo # deviation
Sqo = sum(d**2)/N
Sqi = Sqo/m**2
Sm = N*Sqo/(N*sum(qi**2) - sum(qi)**2)
Sb = Sqo*sum(qi**2)/(N*sum(qi**2) - sum(qi)**2)
return Sm,Sb,Sqi,Sqo
qi = r_[xi,xd]
qo = r_[yi,yd]
Sm,Sb,Sqi,Sqo = regression_analysis(qi,qo)
print('Sm, Sb, Sqi, Sqo = ')
print('%6.4f, %6.4f %6.4f %6.4f' % (Sm,Sb,Sqi,Sqo))
Estimate: y = ax + b N = 22 Degrees of freedom $\nu$ = 20 a = 1.08 $\pm$ 0.006 b = -0.85 $\pm$ 0.037 R^2 = 0.997 Syx = 0.208 y = 1.08 x + -0.85 $\pm$ 0.093 V Sm, Sb, Sqi, Sqo = 0.0002, 0.0063 0.0336 0.0394
Taking \(3 \sigma\) for \(99.7\%\) certainty level:
\(m = 1.08 \pm 0.0006\)
\(b = -0.85 \pm 0.018\)
When we use calibration curve for the measurement:
\(q_i = (q_o + 0.85)/1.08 = 0.9233 q_0 + 0.785\)
\(S_{q_i}^2 = 0.034\)
If we measured for instance 4.32 kPa, we transfer it to 4.79 kPa and write:
(4.32*0.85)/1.08 = 4.787
The final result is¶
\(P = 4.79 \pm 0.102\) kPa (99.7\%)