BookmarkSubscribeRSS Feed
☑ This topic is solved. Need further help from the community? Please sign in and ask a new question.
SimoneStefano96
Fluorite | Level 6

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 

1 ACCEPTED SOLUTION

Accepted Solutions
Rick_SAS
SAS Super FREQ

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;

 

View solution in original post

3 REPLIES 3
sbxkoenk
SAS Super FREQ

Hello,

 

I can only say that the zero-probabilities are not allowed indeed.

 

From the doc:

sbxkoenk_1-1640626166159.png

Thanks,

Koen

Rick_SAS
SAS Super FREQ

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;

 

SimoneStefano96
Fluorite | Level 6

Dear Rick,

 

I am so thankful for your reply, which definetely made my day.

 

Yours sincerely,

 

Simone

sas-innovate-2024.png

 

Time is running out to save with the early bird rate. Register by Friday, March 1 for just $695 - $100 off the standard rate.

 

Check out the agenda and get ready for a jam-packed event featuring workshops, super demos, breakout sessions, roundtables, inspiring keynotes and incredible networking events. 

 

Register now!

Multiple Linear Regression in SAS

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.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 3 replies
  • 416 views
  • 3 likes
  • 3 in conversation