new;cls;
p=0.4;
n=100;
r=100; @100行@
c=3; @生成3列100行二项式分布随机变量@
w=reshape(sumc((rndu(r*c,1).>(0|pbinom(seqa(0,1,n),p,n))')')-1,r,c);
print w;
proc pbinom(x,p,n); @累积概率分布@
local xx,d;
      xx = trunc(x).*(x.>=0); @ trunc取整或者删除小数@
      if maxc(maxc(xx)) ge 0;
         d = cumsumc(dbinom(seqa(0,1,maxc(maxc(xx))+1),p,n));
         d = reshape(d[xx+1],rows(x),cols(x));
      else;
      d = zeros(rows(x),cols(x));
      endif;
retp(d.*(x.>=0));
endp;
proc dbinom(x,p,n); @概率密度@
local xx,d;
    if p==0; d = (x.==0);
    elseif p==1; d = (x.==n);
    else;
    xx = x.*(x.==trunc(x)).*(x.>=0).*(x.<=n);
    d = exp(lnfact(n)-lnfact(xx)-lnfact(n-xx)+xx*ln(p)+(n-xx)*ln(1-p));
  endif; @ lnfact计算阶乘函数的自然对数@
retp(d.*(x.==trunc(x)).*(x.>=0).*(x.<=n));
endp;