function [ b,stdb ] = linregr( X,Y )
% multilinear regression
% Y    =  X  * b     + E
% (Jx1)  (JxP)   (Px1)   (Jx1)
%
% b      (Px1) = fitting coefficients
% E      (Jx1) = absolute error in each fitting point
% stdb   (Px1) = standard deviation of fitting coefficients
% Refs: https://se.mathworks.com/help/stats/coefficient-standard-errors-and-confidence-intervals.html
% or Press et. al. pp 558 ->
J = length(Y);
P = size(X,2);
b = X\Y;        % Least squares solution to the equation Y=X*b
E = Y - X*b;    % Errors in each point
MSE = (E'*E)/(J-P);         % Mean squared errors 
Covar = MSE.*inv(X'*X);     % Covariance matrix
stdb = diag(sqrt(Covar));   % Standard error of fitting coefficients
end

