function P = randbinsym(n)

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

Aeq = is_stat_sym(n);
Aeq = [Aeq;oo'];
Beq = zeros(size(Aeq,1),1); Beq(end)=1;

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

return


function Aeq = is_stat_sym(n)
Aeq = zeros(0,2^n);
rind = 0;
for i=1:n
  rind = rind + 1;
  for xi = 1:2^n
    x = i2c(xi,n);
    
    Aeq(rind,xi) = x(i);
  end
end
return



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