The usual way to think about the problem is to consider the
probability $\overline{p}$ of the opposite event we're interested in:
that there is no birthday"collision" at all within a group of
$k$ people. This, we're said (without really explaining why), is
easier to calculate, and it follows from it that we can retrieve the
probability $p$ of our goal event (i.e. that there is at least
one "collision") as $1 - \overline{p}$.

The above-mentioned essay discusses at length this idea and the
underlying principles, but it also casually challenges the reader with
an interesting side problem: find a way to calculate
$p$ directly, that is, without relying on the probability of
the opposite event. Being a programmer, I decided to tackle this
problem not from the principles, but from the opposite direction, by
trying to derive understanding from "procedural tinkering".

In the first version of this article, I had derived a way of screening the
$k$-permutations
out of the
cartesian power set $N^k$. I thought this was the answer, but someone on
the Cryptography Class discussion forum helped me understand that this
was actually only a rearrangement of the indirect computation (i.e. $p
= 1 - \overline{p}$). The correct way to compute $p$ directly should
rather involve the sum of:

the chance of a collision occurring only on the second choice

the chance of a collision occurring only on the third choice

...

the chance of a collision occurring only on the $k$-th choice

This seemed to make sense, but as I wanted to study this proposition
in more detail on a simplified problem instance, I wrote this Python
program:

which works by enumerating, for every $k$, all the possible trials (of
length $k$) and looking for a "new" collision with each one, "new"
meaning that no subset (of length $< k$) of this particular trial has
ever been causing a collision. The probabilities for each $k$ are then
summed to yield the overall, direct probability of at least one
collision. Note that this code relies on the itertools
module's product
function to generate the cartesian powers of $k$ (i.e. all possible
trials of length $k$) and the
set data structure for easy past and current collision
detection. Once fully convinced that this was working, the obvious
next step was to derive an analytic formulation for it. By studying
some actual values from my program, I figured out that the probability
for a collision occurring only on the $k$-th choice should be:
\[ P_{only}(n, k) = \frac{\left(\frac{n! \cdot (k-1)}{(n-k+1)!}\right)}{n^k}\]
meaning that the total, direct probability of at least one collision
is the sum:
\[ P_{any}(n, k) = \sum_k{\frac{\left(\frac{n! \cdot (k-1)}{(n-k+1)!}\right)}{n^k}}\]
As this wasn't still fully satisfying, because it doesn't yield an
intuitive understanding of what's happening, the same helpful person
on the Cryptography forum offered an equivalent, but much better
rewrite:
\[ P_{only}(n, k) = \left(\frac{n}{n} \cdot \frac{n-1}{n} \cdot \cdot \cdot \frac{n-k+2}{n}\right) \cdot \left(\frac{k-1}{n}\right)\]
which can be easily understood by imagining a bucket of $n$ elements:
the chance of a collision happening exactly at choice $k$ is the probability that
there was no collision with the first $k-1$ choices
($\frac{n}{n} \cdot \frac{n-1}{n} \cdot \cdot \cdot$, increasing as
the bucket fills up) multiplied by the probability of a collision at
$k$, which will happen $k-1$ times out of $n$, if we assume that the
bucket is already filled with $k-1$ different values (by virtue of the
no-previous-collision assumption).

Although the conclusion of my previous attempt was that it's often
possible to derive mathematical understanding from "procedural
tinkering" (i.e. from programming to math), I'm not so sure anymore,
as this second version is definitely a counterexample.