# Random number generation with specific distribution

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 figured with small benchmarking that one problem is the generation of a random number that is needed in breeding. This breeding happens tens of thousands of times, but I tested my method randomNumber thus:

Timing[Do[randomNumber, {128}]]
{1.263, Null}

Timing[Do[randomNumber, {256}]]
{2.199, Null}


This is intolerable.

So, I ask, what distribution I should use and how? I would like it to be 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.

At the moment I use this monstrosity:

randomNumber[n_] := (Floor[InverseFunction[HarmonicNumber, 1, 2][RandomReal[] HarmonicNumber[n, 1/3], 1/3]] // N) + 1


This has all the properties described above, I can even control the slope by changing the numbers in HarmonicNumber[n, 1/3], 1/3]].

But this method is too slow. I am clearly thinking too complex here. Any suggestions?

• Yeah, wasn't sure if this is correct forum. If you want to suggest some other, would happily go there... Jan 21 '13 at 23:58
• Once you find your distribution you can use RandomVariate
– ssch
Jan 21 '13 at 23:59
• Hmmm, maybe I should ask about the distribution from math forum...I can then ask how to implement it quickly here. Jan 22 '13 at 0:02
• Note that the problem you describe (not of the generation of random number, but really the general problem) sounds a lot like a vehicle routing problem. Jan 23 '13 at 12:20
• @JacobAkkerboom You are correct, this is a variant of that problem. I think I will also be doing a project about solving this at my heuristics class. There are also restrictions, like the travellers need to spend varying times at places and each place has a time window where they should arrive. Jan 23 '13 at 16:49

I think you can do this with Quantile. Being unfamiliar with that I'll show something I think is equivalent using Interpolation. The idea is to give an "inverse" based on your weights, and create a lookup function for random reals between 0 and 1.

For example, suppose you want to find random integers from 1 to 10, weighted by inverse of reversing those integers.

len = 10;
weights = Reverse[Range[len]]/Total[Range[len]];
func = Interpolation[
Transpose[{Accumulate@Prepend[weights, 0], Range[0, len]}],
InterpolationOrder -> 0];


Now define a function to compute n random integers according to the desired weights.

randomNumber[n_] := func[RandomReal[{0, 1}, n]]


Bad style, I know, using global values like that. Easy to clean if you like.

Here is an example to give an impression that the weights are honored.

randomNumber

(* Out= {2, 3, 3, 7, 2, 3, 5, 2, 1, 5, 1, 3, 6, 1, 1,
4, 4, 9, 3, 1, 1, 1, 1, 7, 1, 3, 4, 7, 1, 6} *)


Speed check:

Timing[randomNumber[10^5];]

(* Out= {0.210000, Null} *)


Another method is to use RandomChoice with weights. It's probably faster still.

randomNumber2[n_] := RandomChoice[weights -> Range[len], n]


Same example.

randomNumber2

(* Out= {4, 2, 1, 5, 10, 3, 6, 1, 1, 4, 1, 3, 7, 2, 3, 6, 6, 6, 2,
3, 1, 1, 4, 6, 2, 2, 5, 7, 5, 8} *)

Timing[randomNumber2[10^5];]

(* Out= {0., Null} *)


Okay, yeah, that was faster.

• Wonderful, will test tonight with that method. Thanks Jan 22 '13 at 13:07