Perron number distribution A Perron number is a real algebraic integer $\lambda$ that is larger than the absolute value of any of its Galois conjugates. The Perron-Frobenius theorem says that any
non-negative integer matrix $M$ such that some power of $M$ is strictly positive has a
unique positive eigenvector whose eigenvalue is a Perron number.  Doug Lind proved the converse: given a Perron number $\lambda$, there exists such a matrix, perhaps in dimension
much higher than the degree of $\lambda$.  Perron numbers come up frequently in many places, especially in dynamical systems.
My question:

What is the limiting distribution of Galois conjugates of Perron numbers $\lambda$ in
  some bounded interval, as the degree goes to infinity?

I'm particularly interested in looking at the limit as the length of the interval goes to
0.  One way to normalize this is to look at the ratio $\lambda^g/\lambda$, as $\lambda^g$
ranges over the Galois conjugates. Let's call these numbers \emph{Perron ratios}.
Note that for any fixed $C > 1$ and integer $d > 0$, there are only
finitely many Perron numbers $\lambda < C$ of degree $< d$, since there is obviously a bound on the discriminant of the minimal polynomial for $\lambda$, so the question is only interesting when a bound goes to infinity. 
In any particular field, the set of algebraic
numbers that are Perron lie in a convex cone in the product of Archimedean places of the field. For any lattice, among lattice points with $x_1 < C$ that are within this cone, the projection along lines through the origin to the plane $x_1 = 1$ tends toward the uniform
distribution, so as $C \rightarrow \infty$, the distribution of Perron
ratios converges to a uniform distribution in the unit disk (with a contribution for each complex place of the field) plus a uniform distribution
in the interval $[-1,1]$ (with a contribution for each real place of the field).
But what happens when $C$ is held bounded and the degree goes to infinity?  This question seems
related to the theory of random matrices, but I don't see any direct translation from
things I've heard.  Choosing a random Perron number seems very different from choosing 
a random nonnegative integer matrix.
I tried some crude experiments, by looking at randomly-chosen polynomials of a fixed degree whose coefficients are integers in some fixed range except for the coefficient of $x^d$
which is $1$, selecting from those the irreducible polynomials whose largest real root is Perron.  This is not the same as selecting a random Perron number of the given degree
in an interval.  I don't know any reasonable way to do the latter except for small enough $d$ and $C$ that one could presumably find them by exhaustive search.
Anyway, here are some samples from what I actually tried.
First, from among the 16,807 fifth degree polynomials with coefficients in the range -3 to 3, there are $3,361$ that define a Perron number.  Here is the plot of the Perron ratios:
alt text http://dl.dropbox.com/u/5390048/PerronPoints5%2C3.jpg
Here are the results of a sample of 20,000 degree 21 polynomials with coefficients
between -5 and 5.  Of this sample, 5,932 defined Perron numbers:
alt text http://dl.dropbox.com/u/5390048/PerronPoints21.jpg
The distribution decidedly does not appear that it will converge toward a uniform distribution on the disk plus a uniform distribution on the interval.  Maybe the artificial bounds on the coefficients cause the higher density ring.

Are there good, natural distributions for selecting random integer polynomials? Is there a
  way to do it without unduly prejudicing the distribution of roots?

To see if it would help separate what's happening,
I tried plotting the Perron ratios restricted to $\lambda$ in subintervals. For
the degree 21 sample, here is the plot of $\lambda$ by rank order:
alt text http://dl.dropbox.com/u/5390048/CDF21.jpg
(If you rescale the $x$ axis to range from $0$ to $1$ and interchange $x$ and $y$ axes,
this becomes the plot of the sample cumulative distribution function of $\lambda$.)
Here are the plots of the Perron ratios restricted to the intervals $1.5 < \lambda  < 2$
and $3 < \lambda < 4$:
alt text http://dl.dropbox.com/u/5390048/PerronPoints21%281.5%2C2%29.jpg
alt text http://dl.dropbox.com/u/5390048/PerronPoints21%283%2C4%29.jpg
The restriction to an interval seems to concentrate the absolute values of Perron ratios even more. The angular distribution looks like it converges to the uniform
distribution on a circle plus point masses at $0$ and $\pi$. 
Is there an explanation for the distribution of radii?  Any guesses for what it is?
 A: You can ask a similar question about Pisot and Salem numbers. Last year together with my project student Charlie Scarr we were looking, in particular, at a possible connection between distribution of the roots inside the unit circle and the Mahler measure of a polynomial. We did not progress too far but Charlie made some interesting observations which can be found in his report http://www.maths.dur.ac.uk/~dma0mb/projects/C_Scarr.pdf
Closer to your question, I think it would be very interesting to check what kind of picture comes out if one restricts to the polynomials with small Mahler measure. A large list of such polynomials is available at Michael Mossinghoff's page on Lehmer's problem. It is in no way a random sample but the ordering by Mahler measure looks quite natural in this context.
A: I've gained some new perspective on this question, based partly on comments and on Hitachi Peach's answer. Instead of editing the original question, I'll write down some more thoughts as a (partial) answer in the hopes that it will inspire someone with different expertise to say more.
First, after Hitachi Peach's comment following his answer, I tried plotting a picture
of all the answers for a couple two of the simplest situations: quadratics and cubics
with a small value of $C$.
Below is a diagram in the coefficient space for quadratic polynomials. The horizontal axis is the coefficient of $x$ and the vertical axis is the constant.
alt text http://dl.dropbox.com/u/5390048/QuadraticSmallPerron.jpg
The unshaded area in the middle are polynomials whose roots are real with maximum absolute
value 5 and minimum absolute value 1; the left half of this area consists of Perron polynomials. The red lines are level curves of the maximum root.
Here is a similar plot for cubic polynomials, this time showing the region in coefficient space where all roots have absolute value less than 2. 
alt text http://dl.dropbox.com/u/5390048/CubicRootsSmaller2.jpg
Among these are 31 Perron polynomials (where the maximum is attained for a positive real root.  Here are their roots, and the normalized roots (divided by the Perron number):
alt text http://dl.dropbox.com/u/5390048/PerronPoints3%282%29.jpg
alt text http://dl.dropbox.com/u/5390048/PerronPoints3%282%29normalized.jpg
After seeing and thinking about these pictures, it became clear that for polynomials
with roots bounded by $C > 1$, as the degree grows large, the volume in coefficient space
grows large quite quickly with degree, and appears to high volume/(area of boundary) ratio. The typical coefficients become large, and
most of the roots seem to change slowly as the coefficients change, so you don't
bump into the boundary too easily. If so, then
to get a random lattice point within this volume, it should work fairly well to first 
find a random point chosen uniformly in coefficient space, and then move to the
nearest lattice point.
With that in mind, I tried to guess the distribution of roots (invariant by complex
conjugation), choose a random sample of $d$ elements chosen independently from this distribution, generate the polynomial with real coefficients having those roots, round
off the coefficients to the nearest integer, and looking at the resulting roots.
To my surprise, many of the roots were not very stable: the nearest integral polynomial
usually ended up with roots fairly far out of bounds, for any parameter values of
several distributions I tried.  (Note: one constraint on the distribution is that the
geometric mean of absolute values must be an integer $\ge 1$. This rules out the
uniform distribution at least for small values of $C$).
After thinking harder about the stability question for roots (as the coefficients are
perturbed), I realized the importance of the interactions of nearby roots. Whenever
there is a double root, the roots move quickly when coefficients are changed --- i.e.,
the ratio of volume in coefficient space to volume in root space is relatively small.
It's as if nearby roots in effect have a repulsive force. The joint distribution of roots is important: you get wrong answers if you treat them as independent.
With this in mind, I tried an experiment where I clicked on a diagram to put in roots
for a controlling real polynomial by hand, while the computer found the roots of the nearest polynomial
with integer coefficients. With a little practice, this worked well. New roots "prefer" to
be inserted where the existing polynomial is high, so I shaded the diagram by absolute value of the polynomial, to indicate
good places for a new root.  Sometimes, roots of the controlling polynomial become disassociated from roots of the nearest integer polynomial, and the result is often an out-of-bounds root not near any controlling root. In that case, deleting control roots that are disassociated brings it back into line.  As the control roots are moved around, the algebraic integers jump in discrete steps, and these steps are much smaller when the control root distribution is in a good region of the parameter space.
Here's a screen shot from the experiment, (which is fun!):
alt text http://dl.dropbox.com/u/5390048/ControlRoots.jpg
Here, the convention is that each control point above the real axis is duplicated with its complex conjugate. 
Each control point below the real axis is projected to the real axis, and gives a real
root for the control polynomial. All the control roots are shown in black, and the roots of the nearest integer polynomial are shown in red. For these positions, the red roots are nicely associated with black roots.  It is a Perron polynomial, because a
real root has been dragged so that it has maximum modulus.
In the next screenshot, I have dragged several control roots into a cluster around
 11 o'clock. The red roots weren't happy there, so they disassociated from the control roots
and scattered off in different directions, one of them out to a much larger radius. This is a good indication that 
the ratio of coefficient-space volume to root-space volume is small. This polynomial is not
Perron.
alt text http://dl.dropbox.com/u/5390048/RootPerturb-disassociated.jpg
This experiment is much trickier for $C$ close to $1$: the coefficients are much smaller for a given degree, which makes the roots much less stable. They become more stable when there
are lots of roots spread out fairly evenly, mostly near the outer boundary.
Here is one method that in principle should give a nearly uniformly-random choice of
a polynomial with roots bounded by $C$, and I think, by approximating with the nearest 
polynomial having integer coefficients, give a nearly uniform choice of algebraic integer
whose conjugates are bounded by $C$:  Start from any polynomial whose roots are bounded by
$C$, for instance, a cyclotomic polynomial.  Choose a random vector in coefficient space,
and follow a $C^1$ curve whose tangent vector evolves by Brownian motion on the unit sphere.
Whenever a root hits the circle of radius $C$, choose a new random direction in which
the root decreases in absolute value (i.e, use diffuse scattering on the surfaces).
The distribution of positions should converge to the uniform distribution in the
given region in coefficient space.
This method should also probably work to find a random polynomial whose roots are inside
any connected open set, and subject to certain geometric limitations (for instance, it can't be
inside the unit circle) a nearly uniformly random algebraic integer of high degree
whose Galois conjugates are inside a given connected open set.
Of course still more interesting than an empirical simulation would be a
good theoretical analysis.
A: There are many contexts in which one can establish that a family of polynomials, or even a single polynomial, has roots which are distributed according to some measure which approximates the uniform measure on the unit circle.
Let $f(x) = a_N x^N + \ldots + a_1 x + a_0$ be a polynomial with real coefficients, and assume that $a_N \ne 0$. Consider the quantity:
$$L_N(f) = \log \left( \sum |a_k| \right) - \frac{1}{2} \log |a_N| - \frac{1}{2} \log |a_0|.$$
A theorem of Hughes and Nikeghbali (The zeros of random polynomials cluster uniformly near the unit circle, Compositio 144 (1998)) says (informally) that if $L_N(f)$ is small compared to $N$, then the roots are distributed uniformly along the unit circle. An ingredient of their result is a theorem of Erdos-Turan (On the distribution of roots of polynomials, Ann. of Math. (2) 51 (1950)) which (under the same hypothesis) shows that the arguments of the roots are distributed uniformly. Using this, we can answer:

Is there an explanation for the distribution of radii? Any guesses for what it is?

If one takes the coefficients to lie in some bounded range, say $[-5,5]$, then $L_N(f)$ has order $O(\log(N))$, which is certainly $o(N)$. Hence the random polynomials generated in this way have roots which are distributing in the unit circle. On the other hand, the polynomials being chosen are those which have a unique largest root of size at most five and usually around $5$, and hence the renormalized roots will appear to cluster around the circle of radius $1/5$. The clustering becomes more pronounced when one restricts to polynomials where the largest root lies in $[5,5-\lambda]$ for smaller $\lambda$ (Hence Thurston's graph with more restricted intervals are more "ring-like." (There is no point mass along the real axis by the theorem of Erdos-Turan.)

Are there good, natural distributions for selecting random integer polynomials? Is there a way to do it without unduly prejudicing the distribution of roots?

Given a large dimensional probability space, a natural way to generate random elements according to the distribution is to use a random walk Metropolis-Hastings algorithm. Concretely, one can start with a random Perron polynomial with largest root $< 5$, and then perturb the coefficients according to some distribution (say a normal distribution with small variance). If the new polynomial is also Perron with largest root $< 5$, choose this polynomial, or else keep the same polynomial. This Markov process should --- under suitable conditions --- generate random polynomials (in the aggregate) according to the required distribution. For example, let's run this algorithm and compare it to the roots of random monic polynomials with coefficients in $[-5,5]$ which have a unique largest root. The roots of polynomials with coefficients in $[-5,5]$ are illustrated in the first graph, and the result of the Metropolis-Hastings algorithm for all Perron polynomials with largest root less than $5$ is given in the second graph:

Notice that the roots of the first graph cluster around the unit circle, as explained above. However, the roots of the second graph are clustering around the boundary. Indeed, this turns out to reflect reality (see the theorem below).

What is the limiting distribution of Galois conjugates of Perron numbers $\lambda$ in some bounded interval, as the degree goes to infinity?

Let's first consider the easier problem of describing "real Perron polynomials" --- that is, monic polynomials with real coefficients which have a unique largest root $\le r$ (necessarily real). The actual Perron numbers form integral lattice points in this space (although not all lattice points, just the ones corresponding to irreducible polynomials; however, the reducible lattice points form a thin subset). If one fixes the degree $N$ and increases the radius $r$, then (because the regions involved are suitable nice) the lattice points are distributed more or less uniformly inside the space (Davenport's Lemma). However, if one fixes $r$ and lets $N \rightarrow \infty$, it is no longer so clear whether the distribution of lattice points can be approximated by the corresponding real region. Studying lattice points in non-convex regions (or even convex ones) often leads to pretty thorny number  theoretic issues, but one can at least hope (and compare with experiment) that the real geometry gives some indication of the truth of the matter. To this end, one has the following:
Theorem: Fix a real positive integer $r > 0$. Let $\Omega_N$ denote the space of monic polynomials with real coefficients which have a real root $\alpha$ such that $r > \alpha > |\sigma \alpha|$ for all conjugates $\sigma \alpha \ne \alpha$. Then, as $N \rightarrow \infty$, the roots of a random polynomial in $\Omega_N$ are distributed uniformly along the circle $|z| = r$.
The proof of this theorem uses the theorem of Hughes and Nikeghbali mentioned above. The point is that one has to obtain estimates of how the quantities $|a_i|$ vary over the space $\Omega_N$, which reduces, in the end, to evaluating and/or estimating various integrals over $\Omega_N$. For example, consider the example above of degree $21$ polynomials with a largest root $< 5$. The model given by the OP choose $a_{21} \in [-5,5]$. It turns out, however, that the expected value of the constant term $a_{21}$ over $\Omega_{21}$ (with $r = 5$) is
$$ \frac{3^2 \cdot 5^{21} \cdot 7 \cdot 13 \cdot 17 \cdot 19}{2^{17} \cdot 11} \sim 8.748 \times 10^{13}.$$
This is pretty big compared to $[-5,5]$! 
Other statistics
There are a few other interesting probabilities one can compute using integrals over spaces like $\Omega_N$ and related spaces. For example, one can consider all monic polynomials of degree $2N$ with the property that their roots all have absolute value at most one. This defines a compact region of (the coefficient space) $\mathbf{R}^N$, and so it has a natural uniform measure. Then the probability that a random such polynomial has no real roots in the interval $[-1,1]$ (equivalently, is positive on $[-1,1]$, equivalently, has no roots on all of $\mathbf{R}$) is equal to
$$\sim  \frac{2C}{\sqrt{2\pi} \cdot (2N)^{3/8}},$$
where the constant $C$ is equal to
$$C = 2^{-1/24} e^{- 3/2 \cdot \zeta'(-1)} =  1.24514 \ldots.$$
Curiously enough, there is a theorem of Dembo, Poonen, Shao, and Zeitouni that a random polynomial (in the more usual sense) whose coefficients are chosen
with (say) identical normal distributions with zero mean is positive in $[-\infty,\infty]$
 with probability $N^{-b + o(1)}$ and positive in $[-1,1]$ with probability $N^{-b/2 + o(1)}$
 for some universal constant~$b/2$, which 
 they estimate be $0.38 \pm  0.015$. On the other hand, 
the exponent occurring above is $3/8 = 0.375$. Is there any direct relationship between
these theorems? For example, does this suggest that $b = 3/4$? (The result of DPSZ actually holds for a very wide range of distributions with the same undetermined $b$, but it does not apply to our very rigid context.)
Further Remarks
There's some further analysis one can make of the space $\Omega_{N}$, or more generally the space of all monic polynomials with real coefficients whose roots have absolute value at most $r$. For example, one can show that a "random" polynomial of degree $N$, subject to the constraint that all its roots have absolute value at most $r$, will be a Perron polynomial (i.e. have a unique largest root) with probability exactly $1/N$ if $N$ is odd and $1/(N+1)$ if $N$ is even. Some of this analysis will appear in my thesis (I will add a link here when it is written!). As my advisor said, "Any question that Thurston asks is probably worth thinking about."
A: (This is more of an extended comment than an answer.)
You speculate whether imposing artificial bounds on the coefficients imposes a bias on the pictures your are producing. There is reason to think that this is possible. Edelman and Kostlan have some nice results on "random" polynomials,  where a possible candidate for "random" is given by taking the coefficients $a_n$ of a polynomial of degree $d$ to be normally distributed with variance
$\binom{d}{n}$. In this case, they show that the expected number of real
roots is $\sqrt{d}$ (see *1,*2); in contrast to the result of Kac mentioned in the
comments.
Suppose one takes "random" polynomials of large degree $d$, all of whose
coefficients are integers in the fixed interval $[-m,m]$. One guess is that as
$d \rightarrow \infty$, the distribution of roots of this polynomial approaches
the uniform measure on the the unit circle. This may even be relatively easy to prove, I haven't thought so much about it beyond some postage stamp heuristics.
(Certainly the distribution of $|z|$ approaches a point measure at $1$;
this follows (essentially) from Proposition 2.1 of *3.)
Suppose one now restricts to irreducible polynomials which have a unique
largest root of size $|\lambda|$. It might not be too much of a stretch
to imagine that the distrubution of the other roots is otherwise unchanged,
and so lie (roughly) on the circle $|z| = 1/\lambda$. I can't tell to what
extent this is in accordance with your diagrams.
It's also not clear to me to what extent your diagrams depend at all on the
Perron property. What happens if one
considers random degree $21$ polynomials with coefficients in $[-5,5]$
and simply normalizes by the absolute value of the largest root - does one obtain
substantially different pictures?
*1: http://www-math.mit.edu/~edelman/homepage/papers/kac.pdf
*2: http://www-math.mit.edu/~edelman/homepage/papers/roots.pdf
*3: http://www.dtc.umn.edu/~odlyzko/doc/arch/polynomial.zeros.pdf
(warning, PDF file is backwards).
