How I generate random numbers with constraints?

How can I generate three numbers randomly with three constraints? for example: I want a,b,c in range of 0 to 40 and I want these three inequalities to be correct: c-a-b<=-10 b+c-a<=10 b-c-a>=-10

[Walter Roberson]

integer or continuous?

[Walter Roberson]

According to email I was sent, the user wants continuous values, not just integer.

[Walter Roberson]

The list of possibilities is not small:

 {a = 10, b = c, c <= 10, 0 < c},
 {a = 40, b = 40, c <= 10, 0 < c},
 {a = 40, b = c+30, 0 < c, c < 10},
 {a = 40, c = 0, 30 <= b, b <= 40},
 {a = b, c = 10, 10 < b, b < 40},
 {a = c, b = 10, 0 < c, c < 10},
 {a = 20-c, b = 10, 0 < c, c < 10},
 {a = -b+10, c = 0, 0 <= b, b <= 10},
 {a = b-10, c = 0, 10 < b, b < 30},
 {a = b-10, c = 0, 30 < b, b < 40},
 {a = b+10, c = 0, 0 < b, b < 10},
 {a = b+10, c = 0, 10 < b, b < 30},
 {a = c+30, b = 40, 0 < c, c < 10},
 {a = 2*c+20, b = c+30, 0 < c, c < 10},
 {b = 10, c = 0, a <= 20, 0 < a},
 {b = 30, c = 0, 20 <= a, a < 40},
 {b = 40, c = 0, 30 <= a, a < 40},
 {a = 40, 0 < c, b < 40, c < 10, c+30 < b},
 {a = b-c+10, 0 < c, 10 < b, b < c+30, c < 10},
 {a = b-c+10, 0 < c, b < 10, c < 10, c < b},
 {a = b+c-10, 0 < c, 10 < b, b < c+30, c < 10},
 {a = b+c-10, 0 < c, b < 40, c < 10, c+30 < b},
 {a = c-b+10, 0 < c, b < 10, c < 10, c < b},
 {b = 10, 0 < c, a < 20-c, c < 10, c < a},
 {b = 40, 0 < c, a < 40, c < 10, c+30 < a},
 {b = c+30, 0 < c, a < 40, c < 10, 2*c+20 < a},
 {c = 0, 0 < b, a < b+10, b < 10, -b+10 < a},
 {c = 0, 10 < b, a < b+10, b < 30, b-10 < a},
 {c = 0, 30 < b, a < 40, b < 40, b-10 < a},
 {0 < c, 10 < b, a < b-c+10, b < c+30, c < 10, b+c-10 < a},
 {0 < c, a < 40, b < 40, c < 10, c+30 < b, b+c-10 < a},
 {0 < c, a < b-c+10, b < 10, c < 10, c < b, c-b+10 < a}

... so probably the simple rejection method is easiest. Generate three values, see if the conditions are satisfied, if not then try again.

[John D'Errico]

Rejection, in this case, won't be terrible, although not terribly efficient either. If I compare the volume of the domain defined by those constraints, compared to the volume of a cube of side length 40 on each side,

V
V =
       3166.7

40^3
ans =
       64000

(1-V/40^3)*100
ans =
       95.052

that suggests you will need to reject 95.052% of the samples you will generate. So, roughly for every 20 samples you generate, throw away 19 of them. Not great, but computers are good at things like this.

Actually, I suppose integer samples would not be that terrible either. Just generate the set of all integer combinations that satisfy the constraints, then sample randomly from that set.

Answer 1

[John D'Errico]

This is NOT that hard to do. Yeah, I know it seems hard as hell.

The first trick is to build a triangulation that describes the volume in question.

vertices = [0           10            0
          5            5            0
         10            0            0
         10           10           10
         30           40            0
         40           30            0
         40           40            0
         40           40           10];

tessellation = [1     2     4     7
   1     4     5     7
   2     3     4     7
   3     4     6     7
   4     5     7     8
   4     6     7     8];

So 6 tetrahedra, based on the vertex set as given.

How did I find that set? I used some software I wrote many years ago, that allows me to slice up a volume in n-dimensions. In this case, I stared with a cube with [0,40] on each axis. The result is what matters however.

A = [-1 -1 1;-1 1 1;-1 1 -1]'
A =
  -1    -1    -1
  -1     1     1
   1     1    -1

vertices*A
ans =
        -10           10           10
        -10            0            0
        -10          -10          -10
        -10           10          -10
        -70           10           10
        -70          -10          -10
        -80            0            0
        -70           10          -10

As you can see, the first column is always less than or equal to -10, the second column always less than or equal to 10, and the third column always greater than or equal to 10.

So effectively, we need to visualize this set of constraints as a quasi-pyramidal volume in R^3. In the end, all you need are the corners of the volume, then delaunayn could have formed the triangulation.

Next, you want to generate random numbers inside that volume, This too is not quite as difficult as it may seem. Well, easier for me to do, since some years ago, I wrote code that does it for a general tessellated domain.

abc = randsc(10,sc3)
abc =
     28.115       32.711       1.9486
     29.906       28.508        2.841
     16.886       10.254       0.5667
     29.667       27.901       7.9517
     15.706       15.787       3.0161
     9.8288       5.0197       2.9441
     20.225       12.401      0.17255
     3.6989        9.926       1.5614
     30.808       28.612       2.3173
     27.678       25.764       4.1306

So, again, we can test this set of random points.

abc*A
ans =
    -58.877       6.5442       2.6469
    -55.573       1.4424      -4.2396
    -26.573      -6.0654      -7.1988
    -49.617       6.1858      -9.7176
    -28.477       3.0969      -2.9353
    -11.904       -1.865      -7.7531
    -32.453      -7.6514      -7.9965
    -12.064       7.7885       4.6657
    -57.103      0.12142      -4.5133
    -49.311       2.2165      -6.0447

As you see, the three columns of that result satisfy the three requirements.

But, since you don't have my code to use, again, not difficult.

The volumes of the 6 tetrahedra are respectively:

Vi =
     583.33
        500
     583.33
        500
        500
        500

Now, just sample randomly from that set, proportionately to the respective volumes.

Sigh. Its probably easiest just to attach the necessary code, rather than give the complete set of instructions on how to sample from a tessellated domain. But I won't do so, because here, there is no need for that course of action.

The pertinent question seems, is this a one-off problem, and you won't need to do it again? How much would I need to teach you to do?

The simple answer is rejection is easily the easiest way to solve it, but not the most efficient.

abc = rand(10000,3)*40;
A = [-1 -1 1;-1 1 1;1 -1 1]';
k = all(abc*A <= [-10 10 10],2);
abc = abc(k,:);

size(abc)
ans =
 455     3

So in a sample of size 10000, 455 of them were retained. The result is nicely uniform over the domain of interest, as rejection methods always work nicely.

plot3(abc(:,1),abc(:,2),abc(:,3),'o')
box on
grid on

If the rejection rate were much higher, then you would need to use the approach I suggest, of dissecting the domain into a set of tetrahedra, then sampling form them. But a 95% rejection rate is not that bad. Just generate more samples than you need. So, if you wanted 1000 samples, than you will need a little over 20000 initial samples in the cube. 25000 should be safe. But what you need to understand is that rejection methods don't ensure how many samples you need to take to achieve a fixed size result. To get that, you nee to use an approach like that which I described initially.

V
V =
     3166.7

1000/(V/40^3)
ans =
      20211