clear
for mc=1:10000
  n = 100;
  k = 5;
  y = randn(n,1);
  x = [ones(n,1),randn(n,k)];
  
  sumse(1) = sum((y - mean(y)).^2);
  wrk(1) = 1;
  for ii=2:(k+1)
    ind = [1,ii];
    b = regress(y,x(:,ind));
    sumse(ii) = sum((y - x(:,ind)*b).^2);
    wrk(ii) = nnz(ind);
  end
  
  aic = n*log(sumse) + 2*wrk + 2*wrk.*(wrk+1)./(n-wrk-1);
  [dum,imin]=min(aic);  

  %bic = n*log(sumse) + log(n)*wrk;
  %[dum,imin] = min(bic);

  %f1 = (sumse(1)-sumse)./(wrk-1);
  %f2 = sumse./(n-wrk);
  
  %f = f1./f2;
  %f(:,1)=NaN;
  %p = 1 - fcdf(f,wrk-1,n-wrk);
  
  %[dum,imin] = min(p);
  %imin(dum>0.05/k)=1;
 
  np(mc) = imin;
end
mean(np==1)
hist(np,[1:(k+1)])    
