bipoiHist:=proc(n,my,kmin,kmax,w)

        //w=1 alle Werte, w=0 nur my, sigma

        local i,bipf,kmi,kma,li;

        begin

         p:=my/n;

         bipf:=stats::binomialPF(n,p):

         kmi:=round(kmin): kma:=round(kmax):

         i:=kmi:li:=[]:

         werte:=[i,bipf(i)] $ i=kmi..kma;

         for i from kmi to kma do

          li:=li.[[i-0.5,0],[i-0.5,bipf(i)],[i+0.5,bipf(i)],[i+0.5,0]];

         end_for;

         hist:=plot::Polygon2d(li,LineColor=[1,0,0]);

         plot(hist);

         if w=1 then

           return(matrix(float([werte])));

         else return(float(["my",n*p]),

                     float(["sigma",sqrt(n*p*(1-p))]))

         end_if;

        end_proc:

bipoiHist(n,3,0,10,1) $ n=5..15

bipoiHist(10^n,3,0,10,1) $ n=1..3

bipoiHist(2000,3,0,10,1)