function P = randgivenstat(m,n,Py)

[Aeq,Beq] = given_stat(m,n,Py);

P = ones(m^n,1); P = P/sum(P);
zz = P*0 + 1e-12; oo = P*0+1;

Aeq = [Aeq;oo'];
Beq(end+1)=1;

f = randn(m^n,1);
P = linprog(f,[],[],Aeq,Beq,zz,oo);

return



function [Aeq,Beq] = given_stat(m,n,Py)
mn = m^n;

Aeq = zeros(0,mn);
Beq = zeros(0,1);
rind = 0;
for i=1:n
  % enforce constraint: marg(X_i) = Py
  for y = 1:m
    rind = rind + 1;
    Beq(rind,1) = Py(y);
    for xi = 1:m^n
      x = i2c(xi,m,n);
      if x(i) == y
        Aeq(rind,xi) = 1;
      end
    end
  end
end

return


function c=i2c(i,m,n)
dims = repmat(m,[1 n]);
c = ind2coord(i,dims);
return