clear all, close all, clc

% Stima di un valore dato N misurazioni affette da rumore gaussiano e
% relativa varianza dell'errore di stima

SignalPower = 1; % può essere +/- 1
NoisePower = 1;

N = 10; % Number of measurements
flagSimula = 1;
M = 1e4; % Simulated Scenarios

NVect = [1 4 16 64]; % To plot the Figure

%% Theory

TErrorVariance = sqrt(NoisePower)/N;

TRatioMeanEstimator = 1;
TRatioStdDevEstimator = sqrt(2/(M-1));

fprintf('Potenza del rumore NP = %g.\n\n', NoisePower);
fprintf('Teoria:\n\nLa varianza dell''errore di stima con N = %d misure è V = %g.\n', N, TErrorVariance);
fprintf('Il rapporto tra stima della varianza dell''errore con la varianza teorica dello stimatore in M = %d scenari ha:\nMedia = %g.\nDeviazione standard = %g\n\n',...
    M, TRatioMeanEstimator, TRatioStdDevEstimator);

%% Simulated Error

if(flagSimula)
    
    Error = zeros(1,M);
    Estimated = Error;

    for i = 1:M

        X = sqrt(SignalPower)*2*(randi(2) - 1.5);
        Noise = sqrt(NoisePower) * randn(1,N);
        Y = X + Noise; % Scalare più matrice fa vettore

        Estimated(i) = mean(Y);

        Error(i) = Estimated(i) - X;

    end

    MCErrorVariance = var(Error);
    MCRatio = MCErrorVariance/TErrorVariance;

    fprintf('Simulato Monte Carlo:\n\nLa varianza dell''errore di stima con N = %d misure in M = %d scenari è V'' = %g.\n', N,M,MCErrorVariance);
    fprintf('Il rapporto tra stima della varianza dell''errore con la varianza teorica è uguale a R = %g\n', MCRatio);

end

%% Figure

FigHandler = figure('units','normalized','position',[0 0 1 1]);

Legenda = char(' '*ones(1+2*numel(NVect), 20));

X = sqrt(SignalPower)*2*(randi(2) - 1.5);

x = X - 4*sqrt(NoisePower):.001:X + 4*sqrt(NoisePower);

StemHandler = stem(X, 0, 'LineWidth', 2, 'MarkerSize', 12);
hold on, grid on
StemHandler.LineStyle = 'none';
StemHandler.Marker = 'x';
Legenda(1, 1:length('True Value')) = 'True Value';

for iN = 1:numel(NVect)
    
    Noise = sqrt(NoisePower) * randn(1,NVect(iN));
    Y = X + Noise; % Scalare più matrice fa vettore
    
    Estimated = mean(Y);
    
    PlotHandler = plot(x, normpdf(x, Estimated, sqrt(NoisePower/NVect(iN))));
    Legenda(2*(iN-1)+2,1:length(sprintf('N = %d', NVect(iN)))) = sprintf('N = %d', NVect(iN));
    
    StemHandler = stem(Estimated, 0, 'LineWidth', 1.5, 'MarkerSize', 10, 'Color', PlotHandler.Color);
    StemHandler.LineStyle = 'none';
    Legenda(2*(iN-1)+3, 1:length(sprintf('Estimated with N = %d', N))) = sprintf('Estimated with N = %d', N);
    
end

xlabel('x')
ylabel('p_X(x)')

title('PDF del valore vero, date le misure')

hleg = legend(Legenda);
hlet.FontSize = 12;

axis([min(x), max(x), 0, 1.05*normpdf(0,0,sqrt(NoisePower/max(NVect)))])