function checkcondexphat(m,n)

P = pnormdim(rand(m^n,1),1);
F = 0*P;
for t=1:n
  xx = randxx(m,t); Xt0=xx{1};
  xx = randxx(m,1); xhat=xx{1};
  [Yhat,K1] = condexphat1(P,F,Xt0,m,n,xhat);
  [Yhat,K0] = condexphat(P,F,Xt0,m,n,xhat);
  er = totdif(K1,K0)
  if totdif(K1,K0) > 1e-6
    keyboard
  end
end
  
return

function [Yhat,K] = condexphat1(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);
for xint = 1:m^n
  x = ind2coord(xint,dims);
  if veq(x(1:t),Xt0)
    join = P(xint);
    cond = join / pXt0;
    K(xint) = K(xint) + cond;
  end
  
  if veq(x(1:t),[Xt1 xhat])
    join = getPX(P,m,n,[Xt1 x(t+1:end)],[1:t-1, t+1:n]);
    cond = join / pXt1;
    K(xint) = K(xint) - cond;
  end
  
end
return



function [Yhat,K] = 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);
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;
  
end
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