function [Yhat,K,TV] = condexphat(P,F,Xt0,m,n,xhat)
% computes Yhat with \hat x_t in place of x_t -- an upper bound

t = length(Xt0);
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);

tt0 = 1:t;
tt1 = 1:t-1;
Xt1 = Xt0(tt1);

pXt0 = getPXt(Xt0,P,m,n) + realmin;
pXt1 = getPXt(Xt1,P,m,n) + realmin;

Yhat = 0;
K = zeros(m^n,1);
TV = 0;
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt); % xtn = \seq{x}{t+1}{n}
  % first term, simple:
  x = [Xt0 xtn];
  xins = coord2ind(x,dims);
  join = P(xins);
  cond = join / pXt0;
  Yhat = Yhat + cond * F(xins);
  K(xins) = K(xins) + cond;
  c1 = cond;
  % second term, trickier:
  join = 0;
  for xt=1:m
    x = [Xt1 xt xtn];
    xins = coord2ind(x,dims);
    join = join + P(xins);
  end
  cond = join / pXt1;
  x = [Xt1 xhat xtn];
  xins = coord2ind(x,dims);
  Yhat = Yhat - cond * F(xins);
  K(xins) = K(xins) - cond;  
  c2 = cond;
  
  TV = TV + abs(c1-c2);
end
TV = TV/2;
return

function p = getPXt(Xt,P,m,n)
t = length(Xt);
if t==0
  p = 1;
  return
end
dimt = repmat([m],[1 n-t]);
dims = repmat([m],[1 n]);
p = 0;
for xint = 1:m^(n-t)
  xtn = ind2coord(xint,dimt);
  x = [Xt xtn];
  xins = coord2ind(x,dims);
  p = p + P(xins);
end
return