Hello,
I'm carrying out a simulation.
I was wondering why this works (I took it from https://blogs.sas.com/content/iml/2013/08/05/simulate-from-multinomial-distribution.html😞
proc iml; call randseed(4321); /* set random number seed */ c = {"black", "brown", "white"}; prob = {0.5 0.2 0.3}; /* probabilities of pulling each color */ X = RandMultinomial(1000, 100, prob); /* 1000 draws of 100 socks with 3 colors */ print X[colname=c];
while this does not:
PROC IML;
call randseed(4321);
c = {"1", "2", "3","4", "5", "6","7", "8", "9","10", "11", "12","13","14"};
prob = {0 0 0.125 0 0.125 0.25 0 0 0.125 0.125 0 0.125 0 0.125 };
X=RandMultinomial(100, 1000, prob);
print X[COLNAME=c];
QUIT;
The log says "ERROR: Matrix X has not been set to a value.".
I assume RandMultinomial can not simulate properly if one of the probabilities is 0. Nevertheless, I can't remove each 0 from the vector, that would be too time-consuming and inconsistent with my goal, I think.
Is there a way to carry out the simulation in SAS/IML keeping the 0's in the vector?
I know I could use the data step-approach reported here
https://blogs.sas.com/content/iml/2016/03/16/simulate-multinomial-sas-data-step.html
but I would like to keep it as simple as possible.
Thanks,
Simone
You can write a short module that strips out the events that have zero probability. You can return Count=0 for the number of times that those categories occurred:
PROC IML;
/* Return count=0 for categories that have prob=0. */
start RandMult(N, NumTrials, prob);
ncols = prod(dimension(prob));
X = j(N, ncols, 0);
idx = loc(prob>0); /* only use categories for which prob > 0 */
posProb = prob[idx];
X[ ,idx] = RandMultinomial(N, NumTrials, posProb);
return X;
finish;
call randseed(4321);
c = "1":"14";
prob = {0 0 0.125 0 0.125 0.25 0 0 0.125 0.125 0 0.125 0 0.125 };
X=RandMult(20, 1000, prob);
print X[COLNAME=c];
QUIT;
Hello,
I can only say that the zero-probabilities are not allowed indeed.
From the doc:
Thanks,
Koen
You can write a short module that strips out the events that have zero probability. You can return Count=0 for the number of times that those categories occurred:
PROC IML;
/* Return count=0 for categories that have prob=0. */
start RandMult(N, NumTrials, prob);
ncols = prod(dimension(prob));
X = j(N, ncols, 0);
idx = loc(prob>0); /* only use categories for which prob > 0 */
posProb = prob[idx];
X[ ,idx] = RandMultinomial(N, NumTrials, posProb);
return X;
finish;
call randseed(4321);
c = "1":"14";
prob = {0 0 0.125 0 0.125 0.25 0 0 0.125 0.125 0 0.125 0 0.125 };
X=RandMult(20, 1000, prob);
print X[COLNAME=c];
QUIT;
Dear Rick,
I am so thankful for your reply, which definetely made my day.
Yours sincerely,
Simone
Registration is open! SAS is returning to Vegas for an AI and analytics experience like no other! Whether you're an executive, manager, end user or SAS partner, SAS Innovate is designed for everyone on your team. Register for just $495 by 12/31/2023.
If you are interested in speaking, there is still time to submit a session idea. More details are posted on the website.
Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.
Find more tutorials on the SAS Users YouTube channel.