> with(plots):with(plottools):with(stats);with(random);

>

> Ncoins:=3:
p:=1/2:
Trials:=2000:

bin:=[seq(0,i=0..Ncoins)]:
for n from 1 to Trials do
ans:=binomiald[Ncoins,p]();
bin:=subsop(ans+1=bin[ans+1]+1,bin);
od:
bin;
display(seq(rectangle([i-3/2,0],[i-1/2,bin[i]/Trials]),i=1..Ncoins+1),
color=green);P0:=%:

>