I'm confused about one part of your analysis. You write that Q^k_l = Q^(k-1)_(l+1) * (1 - 1/(l+1)), but it's not obvious to me why the multipler, (1 - 1/(l+1)), is correct in general (past the second sweep). I went so far as to write a little bit of code to calculate various values of Q^3_l, exactly, from values of Q^2_(l+1), for an 8-element array, and found that it did not hold. (I assume that you were assuming that the initial array permutations were picked uniformly from the set of all permutations). Have I done something wrong or missed something totally obvious?
I think it's right. One way to convince yourself is just to add counters to the implementation of bubble sort, run it for a larger input, and then verify the probabilities (including the conditional probabilities) match what I claim.
You're correct that every permutation of {1, ..., n} is assumed equally likely as an input.
This way of thinking about it might help:
All that "bubble_sweep" does to the input array is to "right-shift" the left-to-right maxima, and "left-shift" the elements between them. The left-shifted elements aren't re-ordered, so they still have the usual probability of being a left-to-right-maxima.
In a bit more detail:
Let L be the sequence of indices, from left-to-right, of left-to-right maxima in the array. Now in s^{th} sweep, let M be the sequence L concatenated with n - s. Assume M has r elements in total now.
I think you could actually implement a (very, very) inefficient "bubble_sweep" that worked along these lines. I think that makes it a bit clearer that, apart from at left-to-right maxima from a previous iteration, we just have subsequences of the original (uniform random) permutation, so the 1 - 1/(l+1) probability holds.
I hope that's clear(er) and I didn't make too many mistakes. The sun is shining and I must go ride my bike now!
Thanks!