function [a, res, s2, V] = reg_lin(X, Y)
  [m,n]=size(X);
  if m~=1 then X=X'; end
  [m,n]=size(Y);
  if m<n then Y=Y'; end
  a = Y\X';
  res = X' - Y*a;
  s2 = sum(res.*res) / (length(X)-length(a));
  V = s2 * inv(Y'*Y);
endfunction