Subversion Repositories wimsdev

!if $wims_read_parm!=slib_header

 !goto proc

!endif

slib_title=Generation of multinomial random data

slib_parms=3\

 1,m\

  ,n\

 ,[p1,p2, .. pj]

slib_author=Sophie Lemaire et Bernadette PERRIN-RIOU

slib_out=random data of m integers following a multinomial law of parameters  n and p=[p_1,...,p_j] (n is an integer and p_i are positive reals of sum <= 1; if the sum of the p_i is less than 1, 1-sum(p_i) is add to the list.

slib_comment= P(X1=k1,X2=k2,..)=n!/(k1!k2!..)*p1^k1*p2^(k2)... In fact, it is preferable to enter only  n-1  p_i because of the tuncature errors : the last one will be  add.

slib_example=3,6,[1/3,1/3]

:proc

!distribute item $wims_read_parm into slib_M, slib_nn

slib_q=!item 3 to -1 of $wims_read_parm

!default slib_M=1

slib_q=!declosing $slib_q

slib_t=!itemcnt $slib_q

slib_s=!sum x for x in $slib_q

!if $slib_s>1 

!exit

!endif

!if $slib_s>0 and $slib_s<1

slib_q=$slib_q, $[1-$slib_s]

!advance slib_t

!endif

slib_mult=

!for slib_r=1 to $slib_M

slib_s=1

slib_N=$slib_nn

slib_V=

!for slib_i=1 to  $[$slib_t-1]

slib_j=!item $slib_i of $slib_q

!readproc slib/stat/binomial  1,$slib_N, $[$slib_j/$slib_s]

slib_V=!append item $slib_out to $slib_V

!distribute item $[$slib_s-$slib_j], $[$slib_N-$slib_out]  into slib_s,slib_N

!next slib_i

slib_V=!append item $slib_N to $slib_V

slib_mult=!append line $slib_V to $slib_mult

!next slib_r

slib_out=!nonempty line  $slib_mult

slib_out=!trim $slib_out