# Finding a distribution for random number generation

I am writing a program for solving the shortest path in travelling salesman problem, with a twist that there are multiple salesmen who partition the cities among themselves, thus creating two part problem, namely partitioning the graph among the salesmen and solving the salesman problem for each partition. I am using mathematica for this, at the moment.

I am using genetic programming, IE. I generate random partitions combined with order of travel and then I breed the solution so that I weight the better solution more, thus allowing all solutions to be included and the diversity to stay good. The actual breeding itself is quite a complex procedure, but that is a question for another forum and another time.

I have alpha version and it works, amazingly, but it is so slow. This would, hopefully, end in commercial use (a start-up firm asked me to look at this), so I need to speed it up before making the program more complex.

I am looking for a good distribution for generating random numbers with the following properties. It should have at least the following properties: discrete, generate numbers between 1 and some larger n, so that $\sum_{i=1}^n P(i) = 1$ and the probabilities to have the property $P(1)>P(2)>P(3)\cdots P(n)$ and I would also like to be able to control how fast the probabilities get smaller and the distribution to be relatively well know, or at least simple to implement.

At the moment I use this kind of method: my distribution is $P(k)=(1/\sum_{i=1}^n 1/(i^{1/v}))/k^{1/v}$. It looks weird but is a solution to a rather simple attempt: I first ask what is the solution $a$ for $\sum_{i=1}^n a/i^{1/v} = 1$ is. This is $1/\sum_{i=1}^n 1/(i^{1/v})$, thus we get the probability by dividing $a$ with $k^{1/v}$. To generate a discrete number I give a random real number between 0 and 1, $m$, and solve what $s$ must be if we want $\sum_{i=1}^s a/i^{1/v} = m$. Finally I take floor of $s$ and add one and get my random number.

This method is too slow, however, and too complex I think. I would like a simpler distribution with easier method for generating random numbers from it that would also have the properties listed above. Any suggestions? Sorry if the terminology is a bit off.

Are you looking for a particular set of functions $P(k)\ (k\leq n)$ with your chosen properties, or how to find a random number according to an already-chosen distribution function? If the latter, then for any $n$ in a reasonable range the best means of generating a random number is likely to be via cumulative table lookup: at the start of your run build a table of values $CumP()$, where $CumP(k) = \sum_{i=1}^k P(i)$; this can of course easily be done in $O(n)$ time by just setting $CumP(0)=0$ and then iterating $CumP(k) = CumP(k-1)+P(k)$ all the way through your table. Once you've done this, choose a random number $0\leq X\leq 1$ uniformly distributed, then (binary) search your $CumP$ table to find the value $k$ with $CumP(k-1)\leq X \leq CumP(k)$; this $k$ will have the correct distribution. This should be fast enough for your needs as long as $n$ is less than a million or so, though you may want to look into so-called aliasing methods to optimize things even a bit more.
As far as families of functions $P$, what you want depends drastically on what sort of falloff you're looking for. I would start with a power-law distribution, $P(k)\approx \alpha^k$ where $\alpha$ is the controllable parameter that says how fast the probability falls off; this means that event $k$ is exactly $\alpha$ times less likely than event $k-1$, for all $k$. Since $\sum_{i=1}^n \alpha^i = \dfrac{\alpha^{n+1}-\alpha}{\alpha-1} = C$, then the precise formula for $P(k)$ to make the probabilities sum to one is $P(k) = \frac{\alpha^k}{C}$, and $C$ can easily be precomputed from $\alpha$ and $n$. If the falloff in this is too steep, then you could also look into a function of the form $P(k) \approx (n-k+1)^s$ for some $s$; if $s=1$ then this is a linear downslope of probabilities such that $k=1$ is $n$ times more likely than $k=n$, whereas if $s$ is larger then the probability falls off more quickly. Here there are explicit formulas for integer $s$, and in particular when $s=1$ then since the sum of the numbers from $1$ to $n$ is $\frac{n(n+1)}{2}$, then in that case $P(k) = \frac{2(n-k+1)}{n(n+1)}$; for other values of $s$ the solution is more complicated, though, and you're probably best off going back to computing the scaling constant $C$ by just summing once over $k=1\ldots n$ at runtime.
• Hmmmm, thanks for answer. I think I will add some more stuff into my original post soon. As for what I am looking for, I think more or less both of your interpretations are needed, though there are other options than this random selection from a distribution for selecting a pair, like selecting uniformly two sets of $m$ candidates and selecting bests within them and breeding them. Also this is really not crucial for their business (I doubt they would have asked me if it was :D ), so I see this all more of a personal project for testing genetic programming. – Valtteri Jan 22 '13 at 9:28