Create a random matrix?

I want to create a random matrix HM (a,b), it satisfies the following conditions:

1. The value of HM is a nonnegative integer from 0 to c.

2. The total value of the elements in a row is less than or equal to d.

you are way that uses loop to made that.

Thanks!

[Pham]

Thanks for everyone's help. After creating the matrix HM, now I want to create 1 vector whose elements are randomly selected from the HM with probability P, this vector satisfies the above constraints. Thank you!

Answer 1

[Roger Stafford]

HM = zeros(a,b);
for k = 1:a
  t = zeros(c,b);
  t(randperm(b*c,d)) = 1;
  HM(k,:) = sum(t,1);
end

Note: This is not a rejection method.

[Pham]

OK. Thank you very much.

Answer 2

[ES]

This should do!

HM=rand(a,b)*d/a;

[Pham]

The value of HM is a nonnegative integer from 0 to c

Answer 3

[Amit]

@Roger: This is a very interesting approach and got me curious. However I think this will not satisfy the random number scenario as requested in the question.

For example, lets take a scenario for a = 5, b = 4, c = 4 and d = 4. In this scenario, the probability of [1 2 0 0] and [1 3 0 0] are same. However, with the approach suggested here, the chances of [1 2 0 0] > [1 3 0 0]. On few test runs, I did not even see a single 3 in the matrix.

The rejection method will atleast guarantee the equal probability of these two scenarios. I also admit that rejection method will be tedious for long matrixes and take a lot of time for the cases where b >> d.

[Image Analyst]

I was actually a bit surprised we didn't get a comment from Roger about how it's not truly random if you have imposed those constraints (as you've found out). Maybe he thought Pham knew that already. I wouldn't have been surprised to see him launch into a discussion of degrees of freedom or something like that. So, you're right - it can't be truly random. Roger can probably tell you why mathematically (and blow both of us away).

[Amit]

I think even with the constrains, it can be random. I think the algorithm implemented here favors few combinations over others. (I am considering random as picking one of the all possible scenario satisfying the constrain)

For example, if we take only 1 row (a = 1), b = 4 and c = 3 and d = 4. Now a random scenario will be picking randomly n location from 4 columns and filling those with integer combination which satisfy the condition. Lets say for n = 2, picking 2 location out of 4 choices will have same probability which can be filled with equal probability with 2 satisfying integers. i.e [1 1 0 0], [1 0 1 0], [1 2 0 0], [1 3 0 0] etc etc all will have same probability.

A rejection method would be something like this can give more randomness:

a = 5;b = 5;c = 3;d = 4;
HM = zeros(a,b);
tmp1 = zeros(1,b);
for i = 1:a
    tmp = d + 1;
   while (tmp > d)
       tmp1 = randi([0 c],1,b);
       tmp = sum(tmp1);
   end
   HM(i,:) = tmp1;
end

This is a very bad idea for the case where b >> d and randi will rarely generate a combination which will satisfy the condition.

[Image Analyst]

It seems to me, though I'm not a mathematician like Roger, that if you have a series of numbers that are random, and if you require their sum to be some number or not greater than some number, that that constraint could be inconsistent with the originally chosen numbers. For example if the series is 3 numbers and the sum must be less than 30 and you have 15 for the first number and 10 for the second number, that means the third number can't be more than 5, so it can't be truly random. It can't be 10 or 20 or some other random number more than 5. You've lost a degree of freedom.

[Amit]

On further thinking, to do the rejection method for b >> d (code below). Here, n spaces are randomly picked and filled. Doing this allows simply lesser rejections.

a = 5;b = 50;c = 3;d = 10;
HM = zeros(a,b);
tmp1 = zeros(1,b);
tmp2 = 0;
for i = 1:a
    tmp = d + 1;
    tmp2 = randi(d);
   while (tmp > d)
       tmp1 = randi([1 c],1,tmp2);
       tmp = sum(tmp1);
   end
   HM(i,randperm(b,tmp2)) = tmp1;
end

[Amit]

@Image Analyst: You are correct for the case where one has to fill n places (in your example 3 places) with c values with a constrain. However the question is different (atleast the way I see it). In your example, the question is having 3 places with total sum less than 30. So get a randomized scenario, one picks 1 such case from all possible valid situations. SThe probability of such a case would be equal for [1 1 1] and [15 10 5] or even [10 10 10].

[Image Analyst]

Well sure you can cherry pick combinations of "random" numbers that work, but the distribution is different. I don't see it as being any different than the frequently asked question where people want Gaussian distributed random numbers between 1 and 5. Sure you can use randn and get a bunch of numbers in that range. But people answering always point out that they aren't truly normally distributed even though they came from the randn() function because you're not allowing anything less than 1 or more than 5. So if you took a trillion numbers, you wouldn't have tails on the distribution that go to -infinity and +infinity - the tails/distribution would be clipped.

[Roger Stafford]

@Amit & Image Analyst: Yes, you are both right. The method I used very definitely does not give equal probabilities to the different solutions. Mea culpa! It favors the "central" values clustered near d/b over those farther away. Using a rejection method which chooses each number uniformly from 0 to c regardless of their sum and then discards those whose sum is too large will avoid that difficulty. However, for large values of b and c, and small d, this unfortunately can lead to some extremely high rejection rates.

On the other hand based on my experience in writing 'randfixedsum', I think avoiding rejection and at the same time having equal probabilities for each successful set of numbers would be a difficult undertaking. The trouble is that for large values of b, c, and d there is an enormous complexity to the space of all successful combinations. To take a simple example, let b = 3, c = 9, and d = 13. Without the d limit there would be an uncomplicated 10 x 10 x 10 cube of possible combinations, but invoking this d limit cuts this cube at a certain hexagonal face and we are faced with the task of determining which of the 1000 combinations lie inside the volume that is at or behind this hexagonal face and selecting each with equal probability. There are in fact 500 such combinations and it is already a slightly messy task to use a random process to select among them uniformly. With large values for b this problem becomes wildly complicated in the resulting discrete versions of the world of b-dimensional polytopes.