Random number generation with specific distribution

Question

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[50], {128}]]
{1.263, Null}

Timing[Do[randomNumber[50], {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?

Answer 1

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[30]

(* Out[47]= {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[48]= {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[30]

(* Out[51]= {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[52]= {0., Null} *)

Okay, yeah, that was faster.