function [r,xy,v,zp,dfe] = weighted_reg_old2(y,w,varz)
% Calculate weighted average using weighted linear least squares
%
% :Model:
%
%   z_i = 1*zpop + noise
%
% :Inputs:
%
%   **Y:**
%        data matrix (nsub x T)
%
%   **w:**
%        weights
%
%   **varz:**
%        variance of data at each time point (nsub x T)
%
% :Outputs:
%
%   **r:**
%        weighted correlation coeff across columns of y
%
%   **xy:**
%        weighted covariance matrix
%
%   **v:**
%        weighted variance estimates for each column of y
%
%   **zp:**
%        weighted mean of each column of y
%

[m,n] = size(y);

% computation stuff
invxwx = inv(X'*W*X);
pxw = invxwx * X' * W;


W = diag(w);                    % Weight matrix

X = repmat(1,m,1);              % Design matrix - 1 column of all ones to calculate average 



zp = px*y;        % weighted population mean



e = y - repmat(zp,m,1);         % residuals



dfe_v = zeros(n,1);             

R = eye(m) - X*px;    % residual inducing matrix



% Calculate effective degrees of freedom



for i=1:n,



    V = diag(varz(:,i));

    dfe_v(i) = (trace(R*V)^2)/trace(R*V*R*V);       % Satherwaite approximation

    

end;



dfe = mean(dfe_v);               % Calculate average df over all time points. 





MSE = (e'*W*e)/dfe;             % Mean square error



xy = inv(X'*W*X)*MSE;           % Covariance matrix for zp;

xy = 0.5*(xy+xy');              % Remove rounding error



v = diag(xy);                   % Variance for zp

r = xy./sqrt(v*v');             % Correlation matrix for zp





return

  
