% Modelo:
%
%   x[n] = theta + r[n], n=1,2,...,N-1
%   donde r[n] es WGN con varianza sigma2 conocida.
%
% Vector de observaciones: x = [x[0], x[1], ..., x[N-1]] 
%
% Se desea estimar theta con el estimador media muestral:
%   theta_e = sum(x)/N 
%
% Se calculan las prestaciones del estimador mediante Nsim simulaciones de
% Montecarlo
%

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

% Obtiene M realizaciones independientes del vector de observaciones y obtiene la estima
% para cada una de ellas
M=1e6;
theta_e=zeros(M,1); % estimas del estimador
for m=1:M % ciclo a las estimas
    x=theta+sqrt(sigma2)*randn(N,1); % vector de observaciones
    theta_e(m)=sum(x)/N; % estima
end

% Prestaciones estimador
sesgo_theta=mean(theta_e)-theta;
varianza_theta=var(theta_e);
MSE_theta=varianza_theta+sesgo_theta^2;

disp(['theta = ', num2str(theta)])
disp(['Sesgo del estimador : ', num2str(sesgo_theta)])
disp(['Varianza del estimador : ', num2str(varianza_theta)])
disp(['MSE del estimador : ', num2str(MSE_theta)])
disp('      ')

% Histograma del estimador
Nb=100;
paso=(max(theta_e)-min(theta_e))/Nb;
edges=min(theta_e):paso:max(theta_e);
Nh = histc(theta_e,edges);
area=paso*M;
figure(1);bar(edges,Nh/area,'histc')
figure(1);title('Histograma del estimador')
