Accuracy and precision of linear systems

Introduction

1. The linear model that is used for generating data and comparing the output of paired measurements. The measurements are done with a sensor that has no bias and has small random error agains anouther sensor with a bias, scaling factor not equal 1 and larger measurement noise. The model is used to determine the errors in measurements and for evaluating the agreement.
2. The sensor model in Simulink that includes linear function as well as the gain, saturation block, bandlimiter and so on. This model allows for simulating different types of sensors.
3. The sensor model in Simulink that includes 3rd order polonomial nonlinear function as well as the gain, saturation block, bandlimiter and so on. This model allows for obtaining linearity error.

Linear model in Matlab

Evaluating error

1. Concordance coefficient slide

Input_variance=0; %or 1

N=100;

sigma_e1=0.001;

sigma_e2=0.01;

beta0=0.01; %0

beta1=1;%1.01; %1

if Input_variance==0

sigma_x=0;

else

sigma_x=0.2;

end

mu_x=1;

x=mu_x+sigma_x*randn(N,1);

e1=sigma_e1*randn(N,1);

e2=sigma_e2*randn(N,1);

y1=x+e1;

y2=beta0+beta1*x+e2;

D=y2-y1;

CCC = f_CCC([y1 y2],0.05);

Ch1 = ['Pearson coeff is: ',num2str(CCC{1, 1}.pearsonCorrCoeff)];

disp(Ch1)

Pearson coeff is: -0.12897

Ch2 = ['Concordance correlation coeff is: ',num2str(CCC{1, 1}.est)];

disp(Ch2)

Concordance correlation coeff is: -0.013396

figure

plot(x,y1,'o')

hold on, plot(x,y2,'*')

xlabel('Input x')

ylabel('Output')

if Input_variance==0

title('Linear model for \sigma_x=0')

else

Ch3 = ['Linear model for \sigma_x=',num2str(sigma_x)];

title(Ch3)

end

legend('y1','y2')

figure

plot(y1,'o')

hold on,

plot([1 length(y1)],[mean(y1) mean(y1)],'b')

plot(y2,'*r')

plot([1 length(y2)],[mean(y2) mean(y2)] ,'r')

xlabel('Measurement reading number')

ylabel('Output')

if Input_variance==0

title('The output of the linear model for \sigma_x=0')

else

Ch3 = ['The output of the linear model for \sigma_x=',num2str(sigma_x)];

title(Ch3)

end

legend('y_1','mean(y_1)', 'y_2', 'mean(y_2)')

annonation_save('',"Fig1.5.jpg", SAVE_FLAG);

figure

plot(D)

xlabel('Measurement reading number')

ylabel('D_i')

if Input_variance==0

title('Diference D=y_2-y_1 for \sigma_x=0')

else

Ch3 = ['Diference D=y_2-y_1 for \sigma_x=',num2str(sigma_x)];

title(Ch3)

end

Evaluating agreement

% Scatter plot

x2 = min(y2):0.001:max(y2);

figure

plot(y1,y2,'o')

hold on

plot(x2,x2)

if Input_variance==0

title('Scatter plot for \sigma_x=0')

else

Ch3 = ['Scatter plot for \sigma_x=',num2str(sigma_x)];

title(Ch3)

title('Scatter plot for \sigma_x=0.01')

end

% Bland Altman plot

method1=y1';

method2=y2';

meanArray = mean([method1;method2]);

diffArray = method1-method2;

meanOfDiffs = mean(diffArray);

stdOfDiffs = std(diffArray);

confRange = [meanOfDiffs + 2.0 * stdOfDiffs, meanOfDiffs - 2.0 * stdOfDiffs];

figure

plot(meanArray,diffArray,'o');

hold on;

line([min(meanArray) max(meanArray)],[confRange(1) confRange(1)],'Color','red','LineStyle','-.');

line([min(meanArray) max(meanArray)],[confRange(2) confRange(2)],'Color','red','LineStyle','-.');

line([min(meanArray) max(meanArray)],[meanOfDiffs meanOfDiffs],'Color','black','LineStyle','-.');

grid; ylabel('Differences'); xlabel('Means');

text(max(method1),1.1*confRange(1),{'\downarrow \mu+2\sigma'})

text(max(method1),0.9*meanOfDiffs,{'\downarrow \mu'})

text(max(method1),0.9*confRange(2),{'\downarrow \mu-2\sigma'})

%xlim([0.97 1.03])

if Input_variance==0

title('Bland Altman plot for \sigma_x=0')

else

Ch3 = ['Bland Altman plot for \sigma_x=',num2str(sigma_x)];

title(Ch3)

end

hold off;

Sensor model

1. Figure 2.4 from chapter 2 was obtained by setting: Input(:,2) = mu_x+0.1*sin(2*pi*t*1);
2. Saturation in slides is modeled by Input(:,2) = mu_x+1.5*sin(2*pi*t*1);

clear_all_but('SAVE_FLAG');

mu_x=1;

fs=10000;

dur=5*fs; % 5 seconds

t=0:1/fs:dur/fs;

Input(:,1) = t;

Input(:,2) = mu_x+0.1*sin(2*pi*t*1); % mu_x+1.5*sin(2*pi*t*1);

model_name = 'SimpleSensorModel';

open_system(model_name);

% Linear model

set_param('SimpleSensorModel/Linear_Nonlinear function','Coefs','[1.01,0]') % Bias is shown as a separate block and therefor last element needs to be zero

set_param('SimpleSensorModel/Bias','Value','0.005')

set_param('SimpleSensorModel/Noise','Cov','[0.0000001]')

set_param('SimpleSensorModel/Gain','Gain','2')

out=sim(model_name);

figure, plot(out.y_sc)

xlim([0.2,2])

ylim([1.5,2.5]) %([-0.5,5.5])

title('Continous sensor output')

L=ylabel('$y_{sc}$ (V)', 'FontSize',14);

set(L,'Interpreter','Latex');

xlabel('Time (s)')

annonation_save('a)',"Fig1.4a.jpg", SAVE_FLAG);

figure, plot(out.y_sd)

xlim([0.2,2])

ylim([1.5,2.5])

title('Discrete sensor output')

L=ylabel('$y_{sd}$ (V)', 'FontSize',14);

set(L,'Interpreter','Latex');

xlabel('Time (s)')

annonation_save('b)',"Fig1.4b.jpg", SAVE_FLAG);

figure, plot(out.y_tc)

xlim([0.2,2])

ylim([0.5,1.5])

title('Transducer output')

L=ylabel('$y_{tc}$', 'FontSize',14);

set(L,'Interpreter','Latex');

xlabel('Time(s)')

Evaluating nonlinearity

model_name = 'SimpleSensorModel';

open_system(model_name);

Input(:,1) = t;

Input(:,2) = t; % input is modeled as an linear function of t

% NonLinear model

% yt=b1*x+b2*xt.^2+b3*xt.^3;

b0=0.1; b1=1; b2=0.08; b3=-0.015;

set_param('SimpleSensorModel/Linear_Nonlinear function','Coefs','[-0.015, 0.08,1,0]') % Bias is shown as a separate block and therefor last element needs to be zero

set_param('SimpleSensorModel/Bias','Value','0.1')

set_param('SimpleSensorModel/Noise','Cov','0')

set_param('SimpleSensorModel/Gain','Gain','1')

out=sim(model_name);

lin_error=PlotNonlinearity(fs, b0, t, Input(:,2), out.y_tc.Data);

annonation_save('',"Fig1.6.jpg", SAVE_FLAG);

fprintf('The terminal non-linearity error is %.1f %%.\n', 100*lin_error(1));

The terminal non-linearity error is 4.8 %.

fprintf('The best fit non-linearity error is %.1f %%.\n', 100*lin_error(2));

The best fit non-linearity error is 4.2 %.

Exersizes