import mir.random.engine;
size_t s = 10000;
double[6] p =[1/6., 1/6., 1/6., 1/6., 1/6., 1/6.]; // probs must add up to one
auto rv = multinomialVar(s, p);
uint[6] x;
rv(rne, x[]);
assert(x[0]+x[1]+x[2]+x[3]+x[4]+x[5] == s);