Consider a stage in this process where exactly $i$ distinct numbers have already been seen $(0 \le i \lt N).$ "Equiprobable" means that on average, out of every $N$ times this stage is reached, in $i$ cases the next number drawn will be among those seen and in the remaining $j=N-i$ cases it will be a new number. Thus, the expected number of draws to see a new number, given $j$ distinct numbers remain to be seen, must be $N/j.$ (This intuitive result is made rigorous by invoking the Geometric distribution: see the Coupon Collector's Problem.)
The expected number of draws to reach $k$ distinct numbers ($k=1, 2, \ldots, N$) is the sum of these values, starting at $j=N$ (no numbers drawn yet) going down to (and including) $j=N-(k-1):$
$$E[\text{number of draws to reach } k]=\sum_{j=N-(k-1)}^N \frac{N}{j} = N(H_N - H_{N-k})$$
where $$H_N = \sum_{j=1}^N \frac{1}{j}$$ is the $N^\text{th}$ harmonic number. (Of course $H_0=0.$)
A special case is $k=N,$ the number of draws expected to collect all $N$ numbers (the coupon-collector problem), equal to $NH_N.$
Here is a plot of the results for a simulation of length 5000. The heights of the bars are the average numbers of turns observed in the simulation. The red curve is the graph of $N(H_N-H_{N-K}).$ You can see how the time needed to observe a new number increases especially sharply at the very end. This is characteristic of the situation for all $N.$
The agreement between the simulation and the theoretical result is excellent. If you wish to explore this further, here is the R
code.
#
# Simulate the process directly by successive sampling -- no shortcuts.
# Implicitly, at step `i+1` all the previous numbers are re-indexed from `1`
# through `i` so that the test of a new number is fast: it must exceed `i`.
# The output is an array of times at which each new number was observed.
#
collect <- function(N) {
cumsum(sapply(1:N-1, function(i) {
count <- 0
repeat{
count <- count+1
if(sample.int(N, 1) > i) break
}
count
}))
}
#
# Harmonic numbers. See https://mathworld.wolfram.com/HarmonicNumber.html
#
H <- function(N) 0.577215664901532861 + digamma(N+1)
#
# Simulation.
#
N <- 30
x <- replicate(5e3, collect(N))
#
# Plotting.
#
plot(rowMeans(x), type="h", lwd=2, ylab="Expectation", xlab=expression(k),
main=paste("Expected Turns for N =", N)) # The results
curve(N * (H(N) - H(N-x)), add=TRUE, col="Red", lwd=2) # Theoretical values