To generate the probability correctly, you need to generate all permutations of the 9 digits.
F := Permute(1..9):
picks := [seq(add(add(10^(i-1)*f[i+j],i=1..3),j=0..8,3),f = F)]:
rng := (min..max)(picks);
rng := 774 .. 2556
# Find unique values
vals := convert(picks,set):
num := numelems(vals);
num := 199
# Determine entry that occurs the most
A := Array(rng):
for val in picks do A[val] := A[val]+1; end do:
imax := max[index](A);
imax := 1764
Further While it is mostly pointless here, a better technique would be to avoid generating all the picks as a list.
T := table(sparse):
for f in F do
val := add(add(10^(i-1)*f[i+j],i=1..3),j=0..8,3);
T[val] := T[val] + 1;
A downside is that it is a little less convenient to find the max entry and associated value, however, because there are multiple indices with the same entry, we need to do this, regardless.
vmax := max(entries(T,'nolist'));
vmax := 3888
for i in [indices(T,'nolist')] do
if T[i] = vmax then print(i) end if
Those agree with the values you listed.
The Statistics package could be used here
X := RandomVariable(EmpiricalDistribution(picks)):
S := Sample(X,1000):