In an indirect way, this related to graphical models and the problem of estimating events in a given model. Your matrix yields a digraph with parallel edges between each pair of nodes $(i,j)$ labelled $P_{ij}$ and $P_{ji} = 1 - P_{ij}$. Now pick any two nodes $i,j$. The probability of the event $T(i) < T(j)$ has to be summed over all DAGs induced by the sampling process in which $T(i)$ is before $T(j)$ in the topological (and total) order.
This is your classic sum-of-products, with an exponential number of terms in the sum. I'd guess that the problem would be #P-hard (i.e intractable to estimate) unless there's some deeper structure that one can assume or exploit (for example, in how graphical model estimation is easy if a related graph has bounded treewidth).
Update: In response to your comment (I ran out of space in the comment field), there are two slightly half-assed things I can suggest:
- You might want to start with junction-tree like methods to get some ideas for what a convergent procedure might look like. While they are different problems, my suspicion is that much of the problem structure is similar.
- On the theory side, even if the problem is intractable, you might be able to get an approximate answer (with guarantees) using similar ideas (or even a reduction) to the method used to approximate the permanent. That's highly nontrivial though. This article reviews some of the literature an