choice = 4

% 1 histogram
% 2 MC sim
% 3 histogram p=N/2
% 4 MC sim p=N/2
% 5 MC sim, alt estimator
switch choice
    case 1
sig2 = 1; 
N = 30; %no of samples 
p = 10; % no of params
H = randn(N,p); 
theta = randn(p,1); 

MC = 1000; 
s2hat = zeros(MC,1);
for k = 1:MC
    w = sqrt(sig2)*randn(N,1); 
    x = H*theta + w; 
    %thhat = (H\x); 
    s2hat(k) = norm(x-H*(H\x),2)^2/N; 
end
figure(1)
hist(s2hat,MC/10)
xlabel('$\hat{\sigma}^2$','Interpreter', 'latex')
title(['N=' int2str(N) ', p=' int2str(p)])
hold on 
a=axis;plot(sig2*[1 1],a(3:4),'-r')
hold off

    case 2
NN = [ 10 20 50 100 200 500 1000];
mse = zeros(1,length(NN)); 
for ii=1:length(NN)
    N=NN(ii);
    
   H = randn(N,p); 
   theta = randn(p,1); 

MC = 1000; 
s2hat = zeros(MC,1);
for k = 1:MC
    w = sqrt(sig2)*randn(N,1); 
    x = H*theta + w; 
    %thhat = (H\x); 
    s2hat(k) = norm(x-H*(H\x),2)^2/(N); 
end 
mse(ii)= norm(s2hat-sig2,2)^2/MC;
end
crb = 2*sig2^2./NN;
figure(1)
loglog(NN,crb,'-',NN,mse,'*-')
legend('CRB', 'MSE')
xlabel('N'), ylabel('mse of $\hat{\sigma}^2$','Interpreter','latex')
title([' p=' int2str(p)])

   case 3
sig2 = 1; 
N = 100; %no of samples 
p = floor(N/2); % no of params
H = randn(N,p); 
theta = randn(p,1); 

MC = 1000; 
s2hat = zeros(MC,1);
for k = 1:MC
    w = sqrt(sig2)*randn(N,1); 

    x = H*theta + w; 
    %thhat = (H\x); 
    s2hat(k) = norm(x-H*(H\x),2)^2/N; 
end
figure(1)
hist(s2hat,MC/10)
xlabel('$\hat{\sigma}^2$','Interpreter', 'latex')
title(['N=' int2str(N) ', p=' int2str(p)])
hold on 
a=axis;plot(sig2*[1 1],a(3:4),'-r')
hold off

    case 4
NN = [2 3 4 5 10 20 50 100 200];
mse = zeros(1,length(NN)); 
for ii=1:length(NN)
    N=NN(ii);
    p=floor(N/2);
    
   H = randn(N,p); 
   theta = randn(p,1); 

MC = 1000; 
s2hat = zeros(MC,1);
for k = 1:MC
    w = sqrt(sig2)*randn(N,1); 
    x = H*theta + w; 
    %thhat = (H\x); 
    s2hat(k) = norm(x-H*(H\x),2)^2/(N); 
end 
mse(ii)= norm(s2hat-sig2,2)^2/MC;
end
crb = 2*sig2^2./NN;
figure(1)
loglog(NN,crb,'-',NN,mse,'*-')
legend('CRB', 'MSE')
xlabel('N'), ylabel('mse of $\hat{\sigma}^2$','Interpreter','latex')
title(['N=' int2str(N) ', p=' int2str(p)])

    case 5
NN = [2 3 4 5 10 20 50 100 200];
mse = zeros(1,length(NN)); 
mse2 = zeros(1,length(NN)); 
for ii=1:length(NN)
    N=NN(ii);
    p=floor(N/2);
    
   H = randn(N,p); 
   theta = randn(p,1); 

MC = 1000; 
s2hat = zeros(MC,1);
s2hat2 = zeros(MC,1);
for k = 1:MC
    w = sqrt(sig2)*randn(N,1); 
    x = H*theta + w; 
    %thhat = (H\x); 
    s2hat(k) = norm(x-H*(H\x),2)^2/(N); 
    s2hat2(k) = norm(x-H*(H\x),2)^2/(N-p+2); %scaling that minim MSE?
end 
mse(ii)= norm(s2hat-sig2,2)^2/MC;
mse2(ii)= norm(s2hat2-sig2,2)^2/MC;
end
crb = 2*sig2^2./NN;
figure(1)
loglog(NN,crb,'-',NN,mse,'*-',NN,mse2,'.-')
legend('CRB', 'MSE-ML','MSE-alt')
xlabel('N'), ylabel('mse of $\hat{\sigma}^2$','Interpreter','latex')
title(['N=' int2str(N) ', p=' int2str(p)])

    otherwise
        error('wrong choice')

end % switch
