MSI Ex1 Matlab code

Dateiinhalt

%% Exercise 4 - Stochastic behavior of estimators
%%

close all
clear all

N = 1000; % number of samples
M = 200; %number of experiments
v_0 = 10; % true value of the voltage
i_0 = 5; % true value of the current
R = v_0/i_0; % true value of the resistance
sigma = 2; % standard deviation used in the generation of the normal noise
n_min = -1; % minimum noise value when uniform noise is considered
n_max = 1; % minimum noise value when uniform noise is considered
epsilon = 0.01; % threshold use to check the estimation error 

%% Normally distributed noise generation
% First we construct a matrix containing the values of the noise for all
% measurements and experiments. Every column represents a different
% experiment, while each row corresponds to a measurement.

n_i = randn(N,M)*sigma; % current measurement noise
n_v = randn(N,M)*sigma; % voltage measurement noise

% Then we define the values of the measurements

i_k = i_0 + n_i; % current measurements
v_k = v_0 + n_v; % voltage measurements


%% Normally distributed noise - SA estimator
% First we calculate the value of the estimated resistance
R_sa = zeros(N,M);
for k=1:N
    R_sa(k,:) = mean(v_k(1:k,:)./i_k(1:k,:),1);
end

% Plot the estimated resistance
figure
plot(R_sa)
title('SA estimator - Gaussian noise')
ylabel('Resistance')
xlabel('Measurement')

% Calculation of the mean value over the experiments
R_sa_mean = mean(R_sa,2);

% Calculation of the standard deviation over the experiments
R_sa_std = std(R_sa,0,2);

%Calculation of the experimental probability P(|R_*(k)-R|>=epsilon)
Pr_sa = sum(abs(R_sa-R)>=epsilon,2)/M;

figure
subplot(2,2,1)
plot(R_sa_mean)
title('Mean of SA est. - G. n.')
xlabel('Measurement')
ylabel('Mean')

subplot(2,2,2)
plot(R_sa_std)
title('Std dev of SA est. - G. n.')
xlabel('Measurement')
ylabel('Std dev')

subplot(2,2,3)
hist(R_sa,50)
title('Distribution of SA est. - G. n.')
ylabel('Realizations')
xlabel('Value')

subplot(2,2,4)
plot(Pr_sa)
title('Prob. error > epsilon for SA est. - G. n.')
xlabel('Measurement')
ylabel('Probabilty')


%% Normally distributed noise - LS estimator
% First we calculate the value of the estimated resistance
R_ls = zeros(N,M);
for k=1:N
    R_ls(k,:) = mean(v_k(1:k,:).*i_k(1:k,:),1)./mean(i_k(1:k,:).*i_k(1:k,:),1);
end

% Plot the estimated resistance
figure
plot(R_ls)
title('LS estimator - Gaussian noise')
ylabel('Resistance')
xlabel('Measurement')

% Calculation of the mean value over the experiments
R_ls_mean = mean(R_ls,2);

% Calculation of the standard deviation over the experiments
R_ls_std = std(R_ls,0,2);

%Calculation of the experimental probability P(|R_*(k)-R|>=epsilon)
Pr_ls = sum(abs(R_ls-R)>=epsilon,2)/M;

figure
subplot(2,2,1)
plot(R_ls_mean)
title('Mean of LS est. - G. n.')
xlabel('Measurement')
ylabel('Mean')

subplot(2,2,2)
plot(R_ls_std)
title('Std dev of LS est. - G. n.')
xlabel('Measurement')
ylabel('Std dev')

subplot(2,2,3)
hist(R_ls,50)
title('Distribution of LS est. - G. n.')
ylabel('Realizations')
xlabel('Value')

subplot(2,2,4)
plot(Pr_ls)
title('Prob. error > epsilon for LS est. - G. n.')
xlabel('Measurement')
ylabel('Probabilty')


%% Normally distributed noise - EV estimator
% First we calculate the value of the estimated resistance
R_ev = zeros(N,M);
for k=1:N
    R_ev(k,:) = mean(v_k(1:k,:),1)./mean(i_k(1:k,:),1);
end


% Plot the estimated resistance
figure
plot(R_ev)
title('EV estimator - Gaussian noise')
ylabel('Resistance')
xlabel('Measurement')

% Calculation of the mean value over the experiments
R_ev_mean = mean(R_ev,2);

% Calculation of the standard deviation over the experiments
R_ev_std = std(R_ev,0,2);

%Calculation of the experimental probability P(|R_*(k)-R|>=epsilon)
Pr_ev = sum(abs(R_ev-R)>=epsilon,2)/M;

figure
subplot(2,2,1)
plot(R_ev_mean)
title('Mean of EV est. - G. n.')
xlabel('Measurement')
ylabel('Mean')

subplot(2,2,2)
plot(R_ev_std)
title('Std dev of EV est. - G. n.')
xlabel('Measurement')

ylabel('Std dev')

subplot(2,2,3)
hist(R_ev,50)
title('Distribution of EV est. - G. n.')
ylabel('Realizations')
xlabel('Value')

subplot(2,2,4)
plot(Pr_ev)
title('Prob. error > epsilon for EV est. - G. n.')
xlabel('Measurement')
ylabel('Probabilty')