Expected number of rounds for this labeling scheme

Question

Task

Here is an interesting math problem:

Let's say that there are \$n\$ indistinguishable unlabeled objects in a bin. For every "round", pull \$k\$ objects randomly out of the bin with equal probability, and apply a label on each object (Nothing happens to objects already having a label). Then put these \$k\$ objects, now all labelled, back into the bin and repeat the process until all objects have a label applied to them. What is the expected number of rounds needed to label all the objects, given \$n\$ and \$k\$?

Your task is to write a program that solves this exact math problem. Turns out, this is quite hard to calculate by hand and is better left for computers to bash out. I will list three methods of solving this problem below. You do not have to follow any of these three methods.

Methods

Method 1

I have asked this question on Math SE recently as I was also stumped on how to do this problem. One answer suggests the following formula: $$\sum_{j=1}^n (-1)^{j+1} \binom{n}{j} \frac{1}{1-\binom{n-j}{k}/\binom{n}{k}}$$ I have implemented this formula in Python: Try it online!

Method 2

Another answer from the same question provides the following formula: $${\bf e_1}^T \left(I_n - \left({\mathcal{1}}_{j\geq i} {\frac{ {{n-i \choose j-i}} {i \choose k-j+i} }{n\choose k}} \right)_{i=0,j=0}^{n-1} \right)^{-1} {\bf 1}$$ where "\$I_n\$ is the \$n \times n\$ identity matrix, \$\bf e_1\$ is the first standard basis vector of \$\mathbb{R}^n\$ and \${\bf1}\$ is the \$n\$-vector of all ones. And in the definition of the other matrix \$\mathcal{1}\$ is the indicator function, i.e. the \$i,j\$-element is \$0\$ if \$j<i\$ and otherwise that quotient of binomials...," quoting from the linked answer.

I have implemented this formula in Mathematica: Try it online!

Method 3

(This is very likely not the golfiest way for any language, but this is what I used to generate the larger test cases.)

Another way is to utilize the recursive nature of this problem, as the results of each round depends on the results of all the previous rounds. The recursive formula is as follows (\$E\$ is the function that actually calculates the expected value):

$$\begin{align}E(n,k)&=f(n,k,n)\\f(n,k,u)&=\begin{cases}0&\text{if }u=0\\\frac{\displaystyle\binom nk+\sum_{i=1}^{\min(u,k)}{\binom ui\binom{n-u}{k-i}}f(n,k,u-i)}{\displaystyle\binom nk-\binom{n-u}k}&\text{otherwise}\end{cases}\end{align}$$

I have implemented this formula in Python: Try it online!

I/O format and other rules

Your program or function should take in two positive integer values \$n,k\$ with \$n\ge k\$ and output the expected value corresponding to those values. The output can be a floating point value with an absolute/relative error under \$10^{-6}\$ from the true value, though the algorithm should theoretically work for arbitrarily large \$n\$ and \$k\$. In practice, it is acceptable if the program is limited by time, memory, or data-type size. This means that although your program may output inaccurate floating point values for certain inputs (as long as the outputs are within the allowed error as mentioned earlier), it should still theoretically be calculating the exact value of the true output.

Because the expected value is always going to be a rational number, you can also output as an exact fraction (does not need to be simplified and can be improper), if your language supports such a data type, or if you want to output a pair of numbers, the numerator and the denominator. If in doubt, always refer back to the default I/O methods.

Your program should also be deterministic, meaning it should always produce the same value when ran. Namely, no "running a gazillion simulations to get the experimental expected value."

Test Cases

n,k -> output
1,1 -> 1
4,2 -> 3.8
6,2 -> 6.928571428571428571428571429
7,4 -> 3.779411764705882352941176471
7,7 -> 1
8,1 -> 21.74285714285714285714285714
8,5 -> 3.458181818181818181818181818
10,3 -> 9.046212999306305594338048699
12,7 -> 4.212148576244657154897357173
20,2 -> 35.30473485789483826080766262
100,7 -> 72.28020367501059456930677280
1000,69 -> 105.2447804224750397369209265
1256,5 -> 1934.893442395300917211652090

---

This is code-golf, so the shortest program in bytes wins!

Answer 1

05AB1E, 13 10 bytes

-3 thanks to @alphalpha

Ýsc¤/<z¹F¥

Attempt This Online!

Uses a variant of the first method. Ports this Jelly answer for the binomial transform.

Ý      calculate [0,1,...,n]
s      swap so k is on top
c      calculate [0Ck, 1Ck, ..., nCk]
¤      push the last item without poping
/      divide by it
<      subtract one
z      and inverse. Luckily for us, 05AB1E have 1/0 = 0
¹F¥    calculate the binomial transform, by calculating the differences n times.

Answer 2

Charcoal, 48 42 bytes

ＮθＮηＩΣＥθ∕Π⁺±⊕…ιθ⁻θ…⁰η×Π…·¹⁻θι⁻Π⁻ι…⁰ηΠ⁻θ…⁰η

Attempt This Online! Link is to verbose version of code. Explanation: Implements a modification of a formula used by several other answers:

$$ \sum_{i=1}^n\frac{-i\choose 1-i+n}{\left(\frac{i-1\choose k}{n\choose k}-1\right)} $$

$$ = \sum_{i=1}^n\frac{{-i\choose n-i+1}{n\choose k}}{{i-1\choose k}-{n\choose k}} $$

$$ = \sum_{i=0}^{n-1}\frac{{-1-i\choose n-i}{n\choose k}}{{i\choose k}-{n\choose k}} $$

$$ = \sum_{i=0}^{n-1}\frac{\prod_{j=i+1}^n(-j)\prod_{j=0}^{k-1}(n-j)}{\prod_{j=1}^{n-i}j\left(\prod_{j=0}^{k-1}(i-j)-\prod_{j=0}^{k-1}(n-j)\right)} $$

Previous 48-byte answer:

ＮθＮηＩΣＥθ∕×Ｘ±¹ιΠ⁻θ⁺…⁰η…⁰⊕ι×Π…·¹⊕ι⁻Π⁻θ…⁰ηΠ⁻⁻θ⊕ι…⁰η

Attempt This Online! Link is to verbose version of code. Explanation: Implements a modification of Method 1:

$$ \sum_{j=1}^n (-1)^{j-1} \binom n j \frac 1 {1 - \binom {n-j} k / \binom n k} $$

$$ = \sum_{j=1}^n \frac {(-1)^{j-1} \binom n j \binom n k} {\binom n k - \binom {n-j} k} $$

$$ = \sum_{j=1}^n \frac {(-1)^{j-1} \frac {n!} {j! (n-j)!} \frac {n!} {(n-k)!}} {\frac {n!} {(n-k)!} - \frac {(n-j)!} {(n-j-k)!}} $$

$$ = \sum_{j=1}^n \frac {(-1)^{j-1} \prod_{i=n-j+1}^n i \prod_{i=n-k+1}^n i} {\prod_{i=1}^j i (\prod_{i=n-k+1}^n i - \prod_{i=n-j-k+1}^{n-j} i)} $$

$$ = \sum_{j=1}^n \frac {(-1)^{j-1} \prod_{i=0}^{j-1} (n-i) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^j i (\prod_{i=0}^{k-1} (n-i) - \prod_{i=0}^{k-1} (n-j-i))} $$

$$ = \sum_{j=0}^{n-1} \frac {(-1)^j \prod_{i=0}^j (n-i) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^{j+1} i (\prod_{i=0}^{k-1} (n-i) - \prod_{i=0}^{k-1} (n-j-i-1))} $$

It's then possible to factor the \$ (-1)^j \$ into the following product and then switch the operands of the subtraction in the denominator to correct the sign:

$$ = \sum_{j=0}^{n-1} \frac {\prod_{i=0}^j (i-n) \prod_{i=0}^{k-1} (n-i)} {\prod_{i=1}^{j+1} i (\prod_{i=0}^{k-1} (n-j-i-1) - \prod_{i=0}^{k-1} (n-i))} $$

However this only reduces the code to 45 bytes:

 ＮθＮηＩΣＥθ∕Π⁺⁻θ…⁰η⁻…⁰⊕ιθ×Π…·¹⊕ι⁻Π⁻⁻θ⊕ι…⁰ηΠ⁻θ…⁰η

Attempt This Online! Link is to verbose version of code.

Substituting \$ i = n - j - 1 \$ then results in the final formula used by the 42-byte solution.

Answer 3

Tried all three methods, with some minor changes to each for golfiness. In order:

Wolfram Language (Mathematica), 52 bytes

b[i=-Range@#,++i+#].(b[-i,#2]/b@##-1)^-1&
b=Binomial

Try it online!

Shortest by a long shot. $$E(n,k)=\sum_{i=1}^n{-i\choose 1-i+n}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

Wolfram Language (Mathematica), 82 bytes

Tr@Last@Inverse@Table[Boole[i==j]-b[#-i,#2-i+j]i~b~j/b@##,{i,#},{j,#}]&
b=Binomial

Try it online!

Reversed both axes of the matrix, summing the last row instead.

Wolfram Language (Mathematica), 85 83 82 bytes

{n,k}Check[1/Tr[m=b[#,i=Range@#]b[n-#,k-i]],0](n~b~k+m.#0/@(#-i))&@n
b=Binomial

Try it online!

Reversed the domain.

Answer 4

Nekomata, 13 bytes

→r$ÇƆ/←ŗ0ɔ$ᵑ∆

Attempt This Online!

A port of @Command Master's 05AB1E answer. Unfortunately, Nekomata doesn't have 1/0=0.

→r$ÇƆ/←ŗ0ɔ$ᵑ∆
→r              [0,1,...,n]
  $Ç            [0Ck, 1Ck, ..., nCk]
    Ɔ/          [0Ck/nCk, 1Ck/nCk, ..., (n-1)Ck/nCk]
      ←ŗ        decrement and reciprocal
        0ɔ      Append 0
          $ᵑ∆   Take the delta n times

Answer 5

JavaScript (Node.js), 81 bytes

n=>g=(k,j=n,c=(a,b=k)=>!b||c(a,--b)*(a-b)/~b)=>j&&c(n,j)/(c(n-j)/c(n)-1)+g(k,j-1)

Try it online!

Method 1

-5 bytes from tsh

-1 from Arnauld

Answer 6

R, 53 bytes

\(n,k,i=1:n-1)sum((`^`=choose)(-i-1,n-i)/(i^k/n^k-1))

Attempt This Online!

Last few test cases produce incorrect results due to numeric overflow.

Answer 7

Pari/GP, 74 64 bytes

Saved 10 bytes thanks to @alephalpha

---

Use the formula

$$E(n,k)=\sum_{i=1}^n{-i\choose 1-i+n}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

---

Try it online!

b=binomial
E(n,k)=sum(i=1,n,b(-i,1-i+n)*(1/(b(i-1,k)/b(n,k)-1)))

Answer 8

APL (Dyalog APL), 61 bytes

f←{+/(¯1*j+1)×(j!⍺)×÷1-(⍵!⍺-j←⍳⍺)÷⍵!⍺}

Attempt This Online!

not as efficient in bytes as possible, just a translation of formula 1.

Answer 9

Maxima, 71 70 bytes

Saved 1 byte thanks to @att

In Maxima, binomial(-i,1-i+n) can be written in binomial(-i,-1-n) to save byte(s).

---

Use the formula

$$E(n,k)=\sum_{i=1}^n{-i\choose n+1-i}\left(\frac{i-1\choose k}{n\choose k}-1\right)^{-1}.$$

---

Try it online!

E(n,k):=sum(binomial(-i,-1-n)/(binomial(i-1,k)/binomial(n,k)-1),i,1,n)

Answer 10

Jelly, 12 bytes

Ḷc⁹÷c’İ;0I⁸¡

Try it online!

A simple port of @CommandMaster’s 05AB1E answer.