Asymptotics of a Bernoulli-number-like function Tony Lezard asked me the following question which seemed like it should not be too hard but which I did not immediately see how to answer.  Define $f(n,k)$ recursively by $f(1,k) = 1$ and
$$f(n,k) = \frac{1}{1-(1-1/k)^n} \sum_{r=1}^{n-1} \binom{n}{r} \left(\frac{1}{k}\right)^{n-r} \left(1 - \frac{1}{k}\right)^r f(r,k).$$
What can we say about $\lim_{n\to\infty} f(n,k)$?  Even the special case $k=2$ would be interesting:
$$f(n,2) = \frac{1}{2^n-1} \sum_{r=1}^{n-1} \binom{n}{r} f(r,2).$$

Addendum. As explained in Tony Lezard’s answer below, $f(n,k)$ arises as the termination probability of a certain recursive elimination process. The process in the case $n=2$ was studied by Helmut Prodinger, “How to select a loser,” Discrete Math. 120 (1993), 149–159, although Prodinger does not seem to have analyzed $f(n,2)$ specifically.  For other references, see the MO question on The Dance Marathon Problem.
 A: If my computation is correct, then f(n, 2) should be roughly
$$\frac12 \sum_{k \in \mathbb{Z}} 2^{k+s} e^{-2^{k+s}}$$
where s = the fractional part of $\log_2 n$.  (Note the terms of the sum decay rapidly both as $k \to \infty$ and $k \to -\infty$.)

OK, here is the argument.
For n ≥ 2, f(n, 2) is equal to the average value over all subsets S of {1, ..., n} of f(|S|, 2) (including S = {1, ..., n}).  So we may imagine a process where we start with {1, ..., n} and at each step we pick a subset uniformly at random of our current subset.  f(n, 2) is the probability that the last time before we hit the empty set, our set contained just one element.
Consider what happens to each element of {1, ..., n} in this process.  At each stage, if it is still in our subset, we keep it with probability 1/2 and throw it out with probability 1/2.  So, the probability that it remains after r steps is 2-r.  Thus we have n independent variables X1, ..., Xn, each with this exponential distribution, and what we seek is the probability that among them there is a unique maximum.
We may express this probability as a sum over the value of this maximum r.  The probability of this occurring for any given r is
$$n(2^{-r} - 2^{-(r+1)})(1 - 2^{-r})^n = \frac 12 n 2^{-r} (1 - 2^{-r})^n.$$
So, $$f(n,2) = \sum_{r \ge 0} \frac 12 n 2^{-r} (1 - 2^{-r})^n.$$
Write $r = \log_2 n + u$.  Then the rth summand becomes
$$\frac 12 2^{-u} (1 - 2^{-u}/n)^n.$$
From here the derivation is not entirely rigorous, but as $n$ increases, this value is roughly
$$\frac 12 2^{-u} e^{-2^{-u}}.$$
Now u is varying over values of the form $r - \log_2 n$ for $r$ a nonnegative integer, so -u ranges over values of the form $k + s$, s = the fractional part of $\log_2 n$, for $-\infty < k \le \lfloor \log_2 n \rfloor$ an integer.  As n goes to ∞, we obtain the entire sum shown at the top.
A: This is not rigorous, but I think it justifies Michael Lugo's observations that f diverges and the behavior is asymptotically log periodic.
To make the recursion more intuitive, note that f(n,k) is a weighted average of f(1,k)...f(n-1,k). The weight is roughly Gaussian, and it is concentrated near the mode of the binomial expansion of (a + b)^n where a=1/k and b=1-1/k, which happens at about the a^[n/k]b^(n-[n/k]) term, so the terms near f(n-[n/k],k) are heavily weighted. 
The weight drops off rapidly away from the mode m=n-[n/k] as a normal distribution with standard deviation about sqrt(1/k (1-1/k))*sqrt(n). 
Intuitively, the behavior should be similar to a continuous function satisfying 
f(x) = average value of f on the interval [(1-1/k)x - c sqrt(x),(1-1/k)k  + c sqrt(x)].
Fix t large. For x even larger, f(x) can then be expressed as a weighted average (convolution) of f on any interval from [t,t*k/(k-1)]. If you multiply x by k/k-1, then you convolve with a slightly wider function, but it's wider by an exponentially decreasing amount. Thus, the support doesn't spread out to cover all of [t,t*k/(k-1)], and it doesn't become invariant under the remainder when you divide log x by log k/(k-1). 
That just says that some initial conditions would produce oscillations. It doesn't prove that these initial conditions would. 
For t small, the interval [m +- c sqrt(x)] may wrap around [t,t*k/(k-1)], which explains why f almost becomes constant.
A: Douglas, hi, hope you're well. Not a backgammon problem this time but one arising out of a playground game my son was playing recently: n children each select at random one of the k corners of the playground and stand in it. A teacher randomly selects one of the corners and all the children in that corner are eliminated from the game. The game continues in this manner until all children are eliminated or there is a unique winner remaining. f(n,k) is the probability that the game will have a winner. In the actual example k was 4 and n was about 30, but I became interested in the asymptotic behaviour as n grows. I'm surprised that the consensus seems to be that f diverges. Plotting f(n,k) out to about n=1000 shows the log-periodic nature of f but the amplitude diminishes and it looks like it's settling down towards a limit. Hardly rigorous of course!
A: I suspect the limit doesn't exist, believe it or not.  It seems like it should!  But plotting $f(n,2)$ for n from, say, 30 to 2000, it looks like we have
$$ f(n) = C + \omega(n) + o(1) $$
where $C \approx 0.721347$ and $\omega(n)$ is a certain function which is "log-periodic" in the sense that $\omega(n) = \omega(2n)$.  $\omega(n)$ never gets bigger than about $8 \times 10^{-6}$, so you have to look pretty closely to see it.
This sort of very small log-periodic fluctuation appears in certain questions from the mathematics of analysis of algorithms.  For example, the mean length  of the longest run in a random binary word of length $n$ is $\log_2 n + C + P(\log_2 n) + O(n^{-1/2} \log^2 n)$, where $P$ is a continuous periodic function of period 1; see Flajolet and Sedgewick, Analytic Combinatorics, Proposition V.1.  I've also seen similar results for problems involving binary trees.
Something similar seems to be true for $f(n,3)$, except that the log-periodic fluctuations appear to be on the order of $10^{-9}$ and the fluctuations are ``faster''; if we let $\omega$ denote the log-periodic function as above, then we have something like $\omega(n) = \omega(1.5 n)$  (I write 1.5, not $3/2$, because I am not at all confident that the relevant constant is $3/2$; I'm not doing the computations with enough precision to be sure.)
