function [LT,mut]=imputeTSESP(X)
% first row of X are y values and the second row are x values-impute any as
% necessary. 
LT=cov(X'); mut=mean(X'); n=size(X,2); % in our case n=6
ymis=[1 5]; xmis=[4];  % coordinates of first of two missing values, sum of pair is correct.
for j=1:100000
    m=kron(mut',ones(1,n));  % the mean reshaped.
L=(1/n)*(X-m)*(X-m)';   rho2=(LT(1,2)^2)/(LT(1,1)*LT(2,2));
mu=mean(X');
beta=[LT(1,2)/LT(1,1) LT(1,2)/LT(2,2)]; % second coordinate are values of beta_y|x and the first, beta_x|y
for i=ymis
    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)=.5*(t+beta(2)*(X(2,i)-X(2,i+1)));  %use this for EM algorithm
    X(1,i+1)=t-X(1,i);
end
for i=xmis
    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)));
    %X(2,i)=.5*(t+beta(1)*(X(1,i)-X(1,i+1)));  %use this for EM algorithm
    X(2,i+1)=t-X(2,i);
end
%X
LT=(1-1/j)*LT+(1/j)*L;     % this is the cumulative average of the covariance matrices
mut=(1-1/j)*mut+(1/j)*mu;  % cumulative average of the mu values
end