
% Datos:
N=10000;
A=2;
varianza=3;

% Simulacion Montecarlo
Nsim=1e5;
A_e=zeros(Nsim,1);
varianza_e=zeros(Nsim,1);
for ns=1:Nsim
    x=A+sqrt(varianza)*randn(N,1); % Observacion
    A_e(ns)=sum(x)/N;
    varianza_e(ns)=sum((x-A_e(ns)).^2)/N;
%    varianza_e(ns)=sum((x-A_e(ns)).^2)/(N-1);    
end

% Prestaciones estimador de la media
media_A=mean(A_e);
varianza_A=var(A_e);
sesgo_A=media_A-A;
MSE_A=varianza_A+sesgo_A^2;

% Prestaciones estimador de la varianza
media_v=mean(varianza_e);
varianza_v=var(varianza_e);
sesgo_v=media_v-varianza;
MSE_v=varianza_v+sesgo_v^2;

disp(['A = ', num2str(A)])
disp(['Sesgo del estimador de A : ', num2str(sesgo_A)])
disp(['Varianza del estimador de A : ', num2str(varianza_A)])
disp(['MSE del estimador de A : ', num2str(MSE_A)])
disp('      ')
disp(['varianza = ', num2str(varianza)])
disp(['Sesgo del estimador de la varianza : ', num2str(sesgo_v)])
disp(['Varianza del estimador de la varianza : ', num2str(varianza_v)])
disp(['MSE del estimador de la varianza : ', num2str(MSE_v)])

% Histograma A_e
Nb=100;
paso=(max(A_e)-min(A_e))/Nb;
edges=min(A_e):paso:max(A_e);
Nh = histc(A_e,edges);
area=paso*Nsim;
figure(1);bar(edges,Nh/area,'histc')
figure(1);title('Estimador de A')

% Histograma varianza_e
Nb=100;
paso=(max(varianza_e)-min(varianza_e))/Nb;
edges=min(varianza_e):paso:max(varianza_e);
Nh = histc(varianza_e,edges);
area=paso*Nsim;
figure(2);bar(edges,Nh/area,'histc')
figure(2);title('Estimador de la varianza')

