function [yhat , Vygx] = blip(x,mux,muy,Vxx,Vyy,Vyx)
% BLIP(X,MUX,MUY,VXX,VYY,VYX) computes linear predictor
% BLIP inputs:
% x = given data; n by 1 vector
% mux = mean vector for x; n by 1
% muy = mean vector for y; m by 1
% Vxx = covariance matrix for x; n by n
% Vyy = covariance matrix for y; m by m
% Vyx = covariance between y and x; m by n
% BLIP outputs:
% yhat = vector of predicted values for y; m by 1
% Vygx = prediction error covariance matrix; m by m
[n p] = size(x);
[m p] = size(muy);
T = chol(Vxx);
temp = [ x-mux Vyx'];
temp = backsolve(T,temp,1,1);
xi = temp(:,1);
zeta = temp(:,2:(m+1));
yhat = muy + zeta'*xi;
Vygx = Vyy - zeta'*zeta;