The Absent-Minded Passengers Problem: A Motivating Challenge Solved by Computer Algebra

In (S.B. Ekhad and D. Zeilberger, 2020) an exciting case study has been initiated in which experimental mathematics and symbolic computation are utilized to discover new properties concerning the so-called Absent-Minded Passengers Problem. Based on these results, Doron Zeilberger raised some challenging tasks to gain further probabilistic insight. In this note we report on this enterprise. In particular, we demonstrate how the computer algebra packages of RISC can be used to carry out the underlying heavy calculations.


A challenging Email
On January 22, 2020 I received the following email by Doron Zeilberger: Can you (and Sigma) extend theorem 5 of that paper to the general case with k absent-minded passengers?
If you and Sigma can do the fourth moment, and derive the asymptotic in n (with a fixed but arbitrary k), I will donate $100$ to the OEIS in your honor. ...
When I received this email, I was thrilled: First, if one gets such an email -in particular from Doron Zeilberger, one is automatically eager to solve it. And second, I highly appreciate "The On-Line Encyclopedia of Integer Sequences" (OEIS) at http://www.oeis.org and supporting it by a donation of Doron Zeilberger gave an extra strong motivation.
Summarizing, the above email provoked various heavy calculations by means of computer algebra that will be introduced in the following.

The underlying problem and symbolic summation
The combinatorial problem of absent-minded passengers can be introduced as follows.
Definition 2.1. Consider a plane with n ≥ 2 seats and suppose that n passengers enter the plane step-wise taking their seats. In addition, suppose that the first k ≥ 1 passengers are absent-minded, i.e., they lost their seat ticket and take a seat (uniformally) at random. The remaining n − k not absent-minded (but shy) passengers take their dedicated seats (as given in the plane ticket) if it is still free; otherwise, they choose (uniformally) at random one of the still available free seats. For 0 ≤ r ≤ n, let p n,k,r be the probability that exactly r passengers sit in the wrong seat and let X n be the random variable for "the number of passengers sitting in the wrong seat". Then the expected value for the passengers sitting in the wrong seat is E(X n ) = n r=1 r p n,k,r and its variance is The situation of one absent-minded passenger (k = 1) has been considered in [17,5] and has been explored further in [7] for the general case k ≥ 1. Among other fascinating results, closed forms for E(X n ) and V (X n ) have been obtained in [7] by skillful combinatorial arguments. More precisely, the definite sum representations have been derived and simplified to and in terms of the harmonic numbers with the special case H n = S 1 (n)). In the following we will prefer to write such expressions in terms of the modified harmonic numbersS yielding, e.g., the more compact expressions can be simplified to a closed form within seconds by executing the command 1 Remark. Internally, the command EvaluateMultiSums of the package EvaluateMultiSums.m [14] uses the summation paradigms of telescoping, creative telescoping and recurrence solving of the summation package Sigma.m [12]; the underlying algorithms generalize the hypergeometric case [9] to the class of indefinite nested sums and products in the setting of difference fields and rings [8,15,16]. In addition, the calculations are supported by special function algorithms of the package HarmonicSums.m [2,3]. We note that the user is completely dispensed from carrying out the summation tools explicitly. However, all the calculation steps are accomplished by proof certificates (based on the creative telescoping paradigm introduced in [18]). Thus if necessary, a rigorous correctness proof can be extracted.

The generating function approach for large moments
As elaborated in details in [6] the expectation, the variance and higher moments can be elegantly described with the generating function approach. Namely, consider the generating function p n,k,r w r of the probability p n,k,r introduced in Definition 2.1. Then the expectation and variance (see also Definition 2.1) can be straightforwardly connected to the generating function via for l ≥ 1; note that M 1 (n, k) = E(X n ). Then one can define the lth moment by In the above email but also in [6] the task was assigned to compute (besides the known moments m 1 (n, k) = 0 and m 2 (n, k) = V (X n )) further moments m l (n, k) (at least for l = 3, 4).
The combinatorial approach of [7] to derive multi-sum expressions of m l (X n ) for larger l seems hopeless. Even though it would have been pure fun to challenge Sigma.m, similarly as in in [10], with more complicated sums than (2.3). Another exciting approach has been carried out in [6,Theorem 5] by guessing the moments m l (n, k) for the special case k = 1 and l = 2, . . . , 8. As it turns out (and proposed by Doron Zeilberger in his email), this attempt can be pushed forward by the following powerful formula given in [6,Theorem 3]: here (x) k = x(x − 1) . . . (x − k + 1) denotes the Pochhammer symbol. As a consequence, one can calculate straightforwardly any value of M l (n, k) and thus of m l (n, k) for particularly chosen n and k. E.g., we can compute and activating the formula (3.2) for n = 10, k = 2 and l = 3 yields the third moment But even more is possible with the hypergeometric sum representation (3.3). As already promoted in [6] one can derive easily the closed form expressions of M l (X n ) for symbolic n and k. For this task we observe that i.e., for any l ∈ Z ≥0 at most l + 1 summands contribute and the remaining summands vanish with the evaluation w = 1. Moreover the differentiation of the arising building blocks in ( In this way, one can compute M 1 (n, k)(= E(X n )) and rediscovers (2.5). Similarly, one gets, e.g., Note that in the special case k = 1 the factor (−1+k)k (−1+n)n does not contribute; this comes from the fact that in the summation of (3.4) the lower bound k − l should be not negative.
In a nutshell, one can calculate sufficiently many M l (n, k) and using the formula (3.2) one gets the desired moments m l (n, k). E.g., with M 1 (n, k) = E(X n ) given in (2.5) and (3.5) one can reproduce m 2 (n, k) as given in (2.6). Similarly, one obtains, e.g., − 30k 3 + 6k 4 − 24kn + k 2 n + 15k 3 n − 4k 4 n + 16kn 2 − 15k 2 n 2 + 3k 3 n 2 δ(−3 + k) For instance, specializing m 4 (n, k) to k = 1 gives for one absent-minded passenger (as derived in [6,Theorem 5]) and m 4 (n, 2) = 2 − 74 + 13n + 13n 2 (−1 + n)n − 24(1 + 2n)S 2 (n) n S 1 (n) + 48S 2 (n) 2 Following this strategy we succeeded in computing the moments m l (n, k) up to l ≤ 15 on a machine with 12 cores and 1.5TB memory, but failed to proceed due to time and space limitations. Luckily, another great trick from [19] enabled us to compute even more moments. Namely, the calculation of turns out to be much faster (and less memory consuming) than the calculation of (3.4). Given the moments of this so-called exponential moments, one gets back the ordinary moments M l (n, k) by using the following formula from [19]: where S(l, r) denote the Stirling Numbers of the Second kind. In summary, we carried out the following steps: Step 1: Calculate the exponential momentsM l (n, k) with the formula (3.6).
Step 2: Calculate the ordinary moments M l (n, k) with the formula (3.7).
Step 3: Finally, calculate the desired moments m l (n, k) with the formula (3.2). Using this approach, we succeeded in calculating generously 2 the moments m l (n, k) up to l ≤ 21 in a recent amount of time and memory. The moments up to order 16 are online available and the link can be found in [11]. The corresponding timings 3 for the different steps, the total time and the size of the moments are summarized in the table of Figure 1.
l time for step 1 time for step 2 time for step 3 total time size of m l (n, k) 1 0. 16   We note that for l ∈ Z ≥0 the closed form expression m l (n, k) can be given in terms of the harmonic numbers S r (n) and S r (k) (resp.S r (n) andS r (k)) with 1 ≤ r ≤ l. We note further that the sequences generated by the harmonic numbers {S r (n) | r ≥ 1} (resp. {S r (n) | r ≥ 1}) are algebraically independent among the rational sequences, i.e., the representation of m l (n, k) in terms of the harmonic numbers is optimal. Interesting enough, the algebraic independence can be shown with the help of the summation paradigm of parameterized (creative) telescoping; see [13,Example 6.3]. For more general classes of harmonic sums we refer also to [4].

Asymptotic expansions
Two days after Doron Zeilberger's first email, I got a slight update of his challenge: Finally, to get the donation in your honor you have to complete the challenge. Find asymptotic expressions for $m_l(n)$ for at least l=7 (it would be nice to also have l=8), and prove (hopefully automatically) that lim n→∞ m l (n) m 2 (n) l/2 equals 0 if l is odd and (2*l-1)(2*l-3)...1 if l is even. This would be a partial elementary reproof of the asymptotic normality proved for arbitrary r using "fancy" probability of the asymptotic normality, in the Henze-Last paper.
Here we parallelized the calculations on 10 Mathematica kernels. E.g., we needed 53.4 hours to obtain the expansion for m 18 (n, k); however summing up all the timings of the used kernels we needed in total 16 days of CPU time to tackle the case l = 18. The expansions up to order 16 are online available and the corresponding link can be found in [11].
To fully win Doron Zeilberger's challenge, we also dealt with the limit Looking at (4.1) the leading term in the expansion of m 2 (n, k) equals log(n) k. In addition the leading term in m 3 (n, k) is also log(n) k. Thus As a consequence we can confirm that the first values c 2l (k) agree with the double factorial of odd numbers (sequence A001147 in OEIS), i.e., c 2l (k) = (2l − 1)!! = l i=1 (2i − 1) = (2l − 1)! (l − 1)!2 l−1 holds for all 1 ≤ l ≤ 9.

Conclusion
When I received Doron Zeilberger's email, the first calculations were straightforward -except that I was (and I am still not) expert in probability theory and thus made some annoying errors in the beginning. However, pushing the calculations further to the 21 th moment and computing its asymptotic expansions was a challenge. Here the computer algebra packages developed at RISC (Sigma.m and EvaluateMultiSums.m by myself and HarmonicSums.m by Jakob Ablinger) were of great help to calculate these huge expressions. I hope that this little note advertises how existing computer algebra tools in general and in particular those of the combinatorics group at RISC can push forward interesting (combinatorial) research topics. For instance, these packages have been heavily used in the last years to carry out large-scale QCD-calculations in the research field of elementary particle physics; see, e.g., [1] and references therein. Last but not least, I am very grateful to Doron Zeilberger who challenged me with these problems and initiated this fun project.