How to generate a random preference relation on a finite set uniformly?

A preference relation (total preorder, weak order) on a set $S$ is a binary relation $ \precsim $ on $S$ that satisfies:

1. completeness: for any $x, y \in S$, either $x \precsim y$ or $y \precsim x$
2. transitivity: for any $x, y, z \in S$, if $x \precsim y$ and $y\precsim z$, then $x \precsim z$

The number of preference relations on a finite set of cardinality $n$ is given by the ordered Bell number $a(n)$:

$$ a(n) = \sum_{k=0}^n \sum_{j=0}^k (-1)^{k-j} {k \choose j} j^n $$

My question is, how to generate a random preference relation on a finite set uniformly? For example, there are $13$ preference relations on $3$ elements, so each of them should have a probability $1/13$ of being chosen.

Of course, I want the algorithm to have computational complexity as low as possible.

What about a simple recursive algorithm, where you first generate a preference for $n-1$ elements, then put the last one either on one of the equivalence classes, or between two of them, or at first or last position?

-- insert pos elt classList inserts elt
--    between classes pos/2-1 and pos/2 if pos is even, 
--    or in class (pos-1)/2 if pos is odd.
insert 0 elt classes = [elt]:classes
insert 1 elt (firstClass:otherClasses) = (elt:firstClass):otherClasses
insert n elt (firstClass:otherClasses) = firstClass:(insert (n-2) elt otherClasses)


randomPreferences [] = return []
randomPreferences (firstElement:otherElements) = do
  othersPreferences <- randomPreferences otherElements
  position <- range 0 (2 * length othersPreferences - 2)
  return (insert position firstElement othersPreferences)