> m = 10^6; die = 1:6
> s = replicate( m, sum(sample(die, 3, rep=T)) )
> round(table(s)/m, 3)
s
    3     4     5     6     7     8     9    10
0.005 0.014 0.028 0.046 0.069 0.097 0.116 0.125
   11    12    13    14    15    16    17    18
0.125 0.116 0.097 0.069 0.046 0.028 0.014 0.005

> m = 10^6;  pop = 1:365;  n= 25
> match = replicate(m, n - length(unique(sample(pop, n, rep=T))) )
>  mean(match);  mean(match==0)
[1] 0.806556
[1] 0.430406

m = 10^6;  pop = 1:365;  n = 25;  match=numeric(m)
for (i in 1:m) {
    bdays = sample(pop, n, rep=T)
    u = unique(bdays);  nu = length(u)
    match[i] = n - nu  }
mean(match)
## 0.805725
mean(match==0)
## 0.430961