function P = randstatmeas(m,n)

Aeq = is_stat(m,n);

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

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

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


return