%% Ejercicios de estimación paramétrica
% 2022 G.T.A.S. (Grupo Tratamiento Avanzado de Señal). Universidad de Cantabria.
%
%
%% Ejercicio 1:
%
% Modelo: x[n] = A + r[n], \quad n=1,2,...,N-1
% Parametro a estimar: \theta=A
% Estimador: \hat \theta = \frac{a}{N}  \sum_{n=0}^{N-1}x[n]
% Datos:

theta=1; % valor real del parametro a estimar (desconocido)
sigma2=1; % varianza del ruido de medida
N=5; % tamaño de los vectores de observaciones 

a=0.5:0.05:1.0; % Valores del parametro a
Na=length(a);

% Realizaciones de los estimadores:

M=1e6;
theta_e=zeros(M,Na);
for m=1:M % ciclo a las estimas
    x=theta+sqrt(sigma2)*randn(N,1); % vector de observaciones
    for na=1:Na
        theta_e(m,na)=a(na)*sum(x)/N; % estima
    end
end

% Estima sesgo, varianza y MSE para los diferentes valores de a:

sesgo=mean(theta_e)-theta;
varianza=var(theta_e);
MSE=mean((theta_e-theta).^2);

% Representa las prestaciones en funcion de a:

% Representa las prestaciones en funcion de a
figure(10);plot(a,sesgo,'bo-',a,varianza,'ro-',a,MSE,'go-','MarkerSize',10)
figure(10);legend({'sesgo','varianza','MSE'},'FontSize',32,'Location','southeast')
figure(10);grid on
figure(10);xlabel('a', 'FontSize',48)
get(gca, 'XTick');
set(gca, 'FontSize', 24)

% Histograma del estimador con minimo MSE:

[minimo,indice] = min(MSE);
figure(11); histogram(theta_e(:,indice),'Normalization','pdf')
figure(11);title('Histograma del estimador con menor MSE')
figure(11);xlabel('\theta', 'FontSize',48)
get(gca, 'XTick');
set(gca, 'FontSize', 24)

%
%
%% Ejercicio 2:
%
% Modelo: x[n] = s[n] + r[n], \quad n=1,2,...,N-1, donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
% Parametro a estimar: \theta=\sigma^2
% Estimador de máxima verosimilitud: \hat \theta_{ML} = \frac{1}{N}  \sum_{n=0}^{N-1}(x[n]-s[n])^2
%
% Datos:


theta=1; % valor real del parametro a estimar (desconocido)
N=8; % tamanno de los vectores de observaciones 
s=2*ones(N,1); % sennal conocida

% Realizaciones del estimador:

M=1e6;
theta_e=zeros(M,1); 
for m=1:M 
    x=s+sqrt(theta)*randn(N,1); % vector de observaciones
    theta_e(m)=sum((x-s).^2)/N; % estima
end

% Prestaciones del estimador:

sesgo=mean(theta_e)-theta;
varianza=var(theta_e);
MSE=varianza+sesgo^2;

disp(['theta = ', num2str(theta)])
disp(['Sesgo del estimador : ', num2str(sesgo)])
disp(['Varianza del estimador : ', num2str(varianza)])
disp(['MSE del estimador : ', num2str(MSE)])
disp('      ')

% Histograma del estimador:

figure(1); histogram(theta_e,'Normalization','pdf')
figure(1);title('Histograma del estimador')

%
%
%% Ejercicio 3:
%
% Modelo: x[n] = A \ s[n] + r[n], \quad n=1,2,...,N-1, donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
% Parametros a estimar: \theta_1=A, \ \theta_2=\sigma^2
% Estimadores \hat \theta_{1} = \frac{\sum_{n=0}^{N-1}x[n]\ s[n]}{\sum_{n=0}^{N-1} s[n]^2},    \hat \theta_{2} = \frac{1}{N}
% \sum_{n=0}^{N-1}(x[n]- \hat \theta_1 s[n])^2
%
% Datos:

theta1=2; % valor real de la atenuacion
theta2=1;  % valor real de la varianza del ruido
N=20; % tamanno de los vectores de observaciones 
s=ones(N,1); % sennal conocida

% Realizaciones de los estimadores:

M=1e6; 
theta1_e=zeros(M,1); 
theta2_e=zeros(M,1); 
for m=1:M 
    x=theta1*s + sqrt(theta2)*randn(N,1);
    theta1_e(m) = sum(x.*s)/sum(s.^2);
    theta2_e(m )= sum((x-theta1_e(m) * s).^2)/N;
end

% Prestaciones de los estimadores:

sesgo1=mean(theta1_e)-theta1;
varianza1=var(theta1_e);
MSE1=varianza1+sesgo1^2;

sesgo2=mean(theta2_e)-theta2;
varianza2=var(theta2_e);
MSE2=varianza2+sesgo2^2;

disp(['theta1 = ', num2str(theta1)])
disp(['Sesgo del estimador para theta1: ', num2str(sesgo1)])
disp(['Varianza del estimador para theta1: ', num2str(varianza1)])
disp(['MSE del estimador para theta1: ', num2str(MSE1)])
disp('      ')
disp(['theta2 = ', num2str(theta2)])
disp(['Sesgo del estimador para theta2: ', num2str(sesgo2)])
disp(['Varianza del estimador para theta2: ', num2str(varianza2)])
disp(['MSE del estimador para theta2: ', num2str(MSE2)])
disp('      ')

% Histograma de los estimadores:

figure(1); histogram(theta1_e,'Normalization','pdf')
figure(1);title('Histograma del estimador para theta1')
figure(2); histogram(theta2_e,'Normalization','pdf')
figure(2);title('Histograma del estimador para theta2')

%
%% Ejercicio 5:
%
% Modelo: x[n] = A \ s[n] + r[n], \quad n=1,2,...,N-1, donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
% Parametro a estimar: \theta=A
% Estimadores \hat \theta_{ML} =  \hat \theta_{LS} = \frac{\sum_{n=0}^{N-1}x[n]\ s[n]}{\sum_{n=0}^{N-1} s[n]^2}
% Datos:

theta=2; % valor real del parametro a estimar (desconocido)
varr=1; % varianza del ruido
N=20; % tamaño de las observaciones
s=cos(pi/8*(0:N-1)'); % señal conocida

% Realizaciones del estimador:

M=1e6;
theta_e=zeros(M,1); 
for m=1:M 
    x = theta * s + sqrt(varr)*randn(N,1); % vector de observaciones
    theta_e(m) = sum(x.*s) / sum(s.^2); % estima
end

% Prestaciones del estimador:

sesgo=mean(theta_e)-theta;
varianza=var(theta_e);
MSE=varianza+sesgo^2;

disp(['theta = ', num2str(theta)])
disp(['Sesgo del estimador : ', num2str(sesgo)])
disp(['Varianza del estimador : ', num2str(varianza)])
disp(['MSE del estimador : ', num2str(MSE)])
disp('      ')

% Histograma del estimador:

figure(1); histogram(theta_e,'Normalization','pdf')
figure(1);title('Histograma del estimador')

%
%% Ejercicio 6:
%
% Modelo: z[n] = A + s + r[n], \quad n=1,2,...,N-1, donde  r[n] es WGN con varianza \sigma^2.
% Parametros a estimar: \theta_1=s, \quad \theta_2= \sigma^2
% Estimadores \hat \theta_{1} = \frac{1}{N} \sum_{n=0}^{N-1}(x[n]-A), \quad
% \hat \theta_{2} = \frac{1}{N}\sum_{n=0}^{N-1}(x[n]-A- \hat \theta_1)^2
% Datos:

N=20; % longitud de las observaciones
theta1=-0.5; % valor real de theta1 
theta2=3; % valor real de theta2
A = 1;

% Realizaciones de los estimadores:

M=1e5; 
theta1_e=zeros(M,1);
theta2_e=zeros(M,1);
for m=1:M
    x=A+theta1+sqrt(theta2)*randn(N,1); % Observacion
    theta1_e(m)=sum(x-A)/N;
    theta2_e(m)=sum((x-A-theta1_e(m)).^2)/N;
end

% Prestaciones de los estimadores:

sesgo1=mean(theta1_e)-theta1;
varianza1=var(theta1_e);
MSE1=varianza1+sesgo1^2;

sesgo2=mean(theta2_e)-theta2;
varianza2=var(theta2_e);
MSE2=varianza2+sesgo2^2;

disp(['theta1 = ', num2str(theta1)])
disp(['Sesgo del estimador de theta1 : ', num2str(sesgo1)])
disp(['Varianza del estimador de theta1 : ', num2str(varianza1)])
disp(['MSE del estimador de theta1 : ', num2str(MSE1)])
disp('      ')
disp(['theta2 = ', num2str(theta2)])
disp(['Sesgo del estimador de theta2 : ', num2str(sesgo2)])
disp(['Varianza del estimador de theta2 : ', num2str(varianza2)])
disp(['MSE del estimador de theta2 : ', num2str(MSE2)])

% Histogramas de los estimadores:

figure(1); histogram(theta1_e,'Normalization','pdf')
figure(1);title('Histograma del estimador de theta1')
figure(2); histogram(theta2_e,'Normalization','pdf')
figure(2);title('Histograma del estimador de theta2')

%
%
%% Ejercicio 7:
%
% Modelo: x[n] = v\ T \ n + r[n], \quad n=1,2,...,N-1, donde s[n] es una señal conocida y r[n] es WGN con varianza \sigma^2.
% Parametros a estimar: \theta_1=v, \quad \theta_2= \sigma^2
% Estimadores \hat \theta_{1} = \frac{1}{T} \frac{\sum_{n=0}^{N-1}n x[n]}{\sum_{n=0}^{N-1}n^2}, \quad
% \hat \theta_{2} = \frac{1}{N}\sum_{n=0}^{N-1}(x[n]-\hat \theta_1 \ n \ T)^2
% Datos:

N=10; % numero de medidas
theta1=1; % velocidad
theta2=2; % varianza del ruido
T=0.1; % periodo de medida

% Realizaciones de los estimadores:

M=1e5; % numero de realizaciones
theta1_e=zeros(M,1);
theta2_e=zeros(M,1);
n=(1:N)';
for m=1:M
    x = theta1*n*T + sqrt(theta2)*randn(N,1); % Observacion
    theta1_e(m) = (1/T)*sum(x.*n)/sum(n.^2);
    theta2_e(m) = sum((x-theta1_e(m)*n*T).^2) / N;
end

% Prestaciones de los estimadores:

sesgo1=mean(theta1_e)-theta1;
varianza1=var(theta1_e);
MSE1=varianza1+sesgo1^2;

sesgo2=mean(theta2_e)-theta2;
varianza2=var(theta2_e);
MSE2=varianza2+sesgo2^2;

disp(['theta1 = ', num2str(theta1)])
disp(['Sesgo del estimador de theta1 : ', num2str(sesgo1)])
disp(['Varianza del estimador de theta1 : ', num2str(varianza1)])
disp(['MSE del estimador de theta1 : ', num2str(MSE1)])
disp('      ')
disp(['theta2 = ', num2str(theta2)])
disp(['Sesgo del estimador de theta2 : ', num2str(sesgo2)])
disp(['Varianza del estimador de theta2 : ', num2str(varianza2)])
disp(['MSE del estimador de theta2 : ', num2str(MSE2)])

% Histogramas de los estimadores:

figure(1); histogram(theta1_e,'Normalization','pdf')
figure(1);title('Histograma del estimador de theta1')
figure(2); histogram(theta2_e,'Normalization','pdf')
figure(2);title('Histograma del estimador de theta2')

%
%% Ejercicio 9:
%
% Modelo: x[n] = A + r[n], \quad n=1,2,...,N-1, donde r[n] es WGN con varianza diferente en cada muestra: \sigma_n^2.
% Parametro a estimar: \theta=A
% Estimadores \hat \theta_{LS} = \frac{1}{N} \sum_{n=0}^{N-1} x[n], \quad
% \hat \theta_{ML} = \frac{\sum_{n=0}^{N-1} x[n] \sigma_n^{-2}}{\sum_{n=0}^{N-1}\sigma_n^{-2}}
% Datos:

theta=1; % valor real del parametro a estimar
varianzasr=[1 1/3 1/6]'; % varianzas de ruido de los sensores
N=length(varianzasr); % tamaño del vector de observaciones (numero de sensores)

% Realizaciones del los estimadores:

M=1e6;
theta_LS=zeros(M,1);
theta_ML=zeros(M,1);
for m=1:M
    x = theta + sqrt(varianzasr).*randn(N,1); % Observacion
    theta_LS(m) = (1/N)*sum(x);  % estima LS
    theta_ML(m) = sum(x./varianzasr) / sum(1./varianzasr); % estima ML
end

% Prestaciones estimadores:

sesgo_LS=mean(theta_LS)-theta;
varianza_LS=var(theta_LS);
MSE_LS=varianza_LS+sesgo_LS^2;

sesgo_ML=mean(theta_ML)-theta;
varianza_ML=var(theta_ML);
MSE_ML=varianza_ML+sesgo_ML^2;

disp(['theta = ', num2str(theta)])
disp(['Sesgo del estimador LS : ', num2str(sesgo_LS)])
disp(['Sesgo del estimador ML : ', num2str(sesgo_ML)])
disp(['Varianza del estimador LS : ', num2str(varianza_LS)])
disp(['Varianza del estimador ML : ', num2str(varianza_ML)])
disp(['MSE del estimador LS : ', num2str(MSE_LS)])
disp(['MSE del estimador ML : ', num2str(MSE_ML)])

% Histogramas de los estimadores:

figure(1); histogram(theta_LS,'Normalization','pdf')
figure(1);title('Histograma del estimador LS')
figure(2); histogram(theta_ML,'Normalization','pdf')
figure(2);title('Histograma del estimador ML')

%%