function L=imputeTSE&SP(L,X)
% first row of X are y values and the second row are x values
ymis=[1 5]; xmis=[4];  % coordinates of firs of two missing values, sum of pair is correct.
L=cov(X');   rho2=(L(1,2)^2)/(L(1,1)*L(2,2));
beta=[L(1,2)/L(1,1) L(1,2)/L(2,2)]; % second coordinate are values of beta_y|x and the first, beta_x|y
for (i in xmis)
    t=X(1,i)+X(1,i+1);
    X(1,i)=normrnd(.5*(t+beta(2)*(X(2,i)-X(2,i+1))), sqrt(.5*(1-rho2)*L(1,1)));
    X(1,i+1)=t-X(1,i);
end
for (i in ymis)
    t=X(2,i)+X(2,i+1);
    X(2,i)=normrnd(.5*(t+beta(1)*(X(1,i)-X(1,i+1))), sqrt(.5*(1-rho2)*L(2,2)));
end
L=cov(X');