How can I randomly draw an ensemble of unit vectors that sum to zero? Inspired by this question, I would like to determine the probability that a random knot of 6 unit sticks is a trefoil.  This naturally leads to the following question:
Is there a way to sample uniformly from the set of ensembles of $n$ unit vectors $\{ v_i \}_{i=1}^n$ in $\mathbb{R}^d$ that sum to zero?  I would like some sort of analytic expression for the distribution (something I might be able to prove a theorem with), but also an algorithmic process to implement in code.
UPDATE:  It looks like what I want is basically the Hausdorff measure of an algebraic variety.  Can I use the construction of this measure produce an analytic expression for the distribution?
 A: Cantarella and Shonkwiler have now posted a paper on the equilateral random polygon case. Their paper seems to cover a lot of what you're interested in, including a new Markov chain algorithm as well as a new lower bound of 3/4 for the probability that a random equilateral hexagon is un-knotted.
That bound is really a huge underestimate, though -- they report that their Markov chain gives about a probability of $(1.3\pm0.2)\times10^{-4}$ for knottedness.

Older post
Cantarella, Deguchi and Shonkwiler have recently begun investigating various ensembles of random polygons which, while not of the particular form you specified, may still be of interest.  
I can't summarize their approach better than the abstract of their paper:

We build a new probability measure on closed space and plane polygons. The key construction is a map, given by Knutson and Hausmann using the Hopf map on quaternions, from the complex Stiefel manifold of 2-frames in n-space to the space of closed n-gons in 3-space of total length 2. Our probability measure on polygon space is defined by pushing forward Haar measure on the Stiefel manifold by this map. A similar construction yields a probability measure on plane polygons which comes from a real Stiefel manifold.
  The edgelengths of polygons sampled according to our measures obey beta distributions. This makes our polygon measures different from those usually studied, which have Gaussian or fixed edgelengths. One advantage of our measures is that we can explicitly compute expectations and moments for chordlengths and radii of gyration. Another is that direct sampling according to our measures is fast (linear in the number of edges) and easy to code.
Some of our methods will be of independent interest in studying other probability measures on polygon spaces. We define an edge set ensemble (ESE) to be the set of polygons created by rearranging a given set of n edges. A key theorem gives a formula for the average over an ESE of the squared lengths of chords skipping k vertices in terms of k, n, and the edgelengths of the ensemble. This allows one to easily compute expected values of squared chordlengths and radii of gyration for any probability measure on polygon space invariant under rearrangements of edges. 

In a follow-up paper Cantarella, Grosberg, Kusner and Shonkwiler compute the total expected curvature for these random $n$-gons and thus extract bounds on knotting probabilities of hexagons and heptagons.

It turns out there is a relatively large literature on "random equilateral polygons", which is the model you are interested in. I don't think there is any known analytic formula for the probability - you would be looking for the measure of some rather complicated semi-algebraic sets corresponding to the components of the configuration space which correspond to knotted polygons. Thus people have primarily been studying properties of such random knots with Monte Carlo sampling.
Ken Millett is one of the experts on random knots and would be a great person to chat with about these questions. Here's one of his papers titled "The Generation of Random Equilateral Polygons"; which contains a discussion of methods used in the literature following variations of the idea that Joseph O'Rourke described above.
Of course, the hexagonal case has been studied numerically but I couldn't find a precise estimate in my searches. An early paper by Ken Millett which shows it as a data point on a curve of knotting probability versus number of edges can be found here.
A: Here is one approach.  Start with a regular $n$-gon in the $xy$-plane with unit edge lengths;
say its vertices are $v_i$, $i=0,\ldots,n-1$.
Now iterate the following process.
Select a random diagonal, $v_i v_j$. Rotate the chain $v_i, v_{i+1}, \ldots, v_j$
(indices appropriately mod $n$)
as a rigid unit about
the line through $v_i v_j$, by a random angle $\theta \in [0,2\pi)$.
Continue until there is sufficient "mixing."
I illustrate the process below for 30 iterations applied to a hexagon.

          


(Apologies for the scale—the chain wanders away from its initial locaion.)
