function [Mstar,Vstar]=...
    getpost_theta(PSI,X,Ystar,rankRP,sig22,sig12)
Vstar=invpd(kron(sig22,PSI'*X'*X*PSI));
        trm1=PSI'*X'*Ystar(:,1:rankRP)*sig12;
        trm2=PSI'*X'*X*PSI*sig12;
        trm3=PSI'*X'*Ystar(:,rankRP+1:end)*sig22;
        Mstar=Vstar*vec(trm1-trm2+trm3);