2

I'm trying to multiply (element wise) a vector V of length N by a randomly generated number in the range (a,b), while keeping the sum of the vector equal to a total amount, E. I want to do this in MATLAB, but I'm not sure how. Getting random numbers between a certain range I know how to do:

minrand = 0;
maxrand = 1;
randfac = (maxrand-minrand).*rand(1,N) + minrand;

But yeah, beyond that I'm pretty clueless. I guess the random numbers can't really be generated like this, because if we call the random numbers the vector R, then I want that R_1*V1 + R_2*V2 .... + R_N*V_N = E. So I guess it's a big equation. Is there any way to solve it, while putting constraints on the max and min values of R?

user129412
  • 765
  • 1
  • 6
  • 19
  • see this post: http://stackoverflow.com/questions/5622608/choosing-n-numbers-with-fixed-sum – Rashid Dec 06 '14 at 15:43
  • Hm, that does seem related, but that is for a vector to add up to a certain number, not an inner product. I don't see how to generalize it straight away.. An alternative idea I had is perhaps just drawing the random numbers, computing the product, and changing values bigger than (maxrand-minrand)/2 until the sum is close to E, within a certain threshold. This will screw up the uniformity of the random numbers, but its something. – user129412 Dec 06 '14 at 15:54
  • You need to _think_ which distribution you want for the random numbers. Saying they should be _random_ and fulfill the sum condition is not enough – Luis Mendo Dec 06 '14 at 16:23
  • Okay, well for my purposes I'd want them to be as uniformly distributed between minrand and maxrand as possible. In practice this will not happen as they need some fine tuning, but for my application this is ok. – user129412 Dec 06 '14 at 16:24

4 Answers4

2

You can pick pairs of two elements (in all combinations) and add and subtract an equal random number.

% Make up a random vector
N=10;
randfac = 10*rand(1,N);

%OP Answer here:  Given randfac with sum E re-randomize it
E = sum(randfac);
minrand = 0;
maxrand = 2;

disp(randfac)
% v = [6.4685    2.9652    6.6567    1.6153    7.3581    0.0237    7.1025
% 3.2381    1.9176    1.3561]
disp(sum(randfac))
% E = 38.7019

r = minrand + (maxrand-minrand)*rand(N*N,1);

k = 1;
for i=1:N
    for j=1:N
        randfac(i) = randfac(i)-r(k);
        randfac(j) = randfac(j)+r(k);
        k = k + 1;
    end
end

disp(randfac)
% v = [5.4905    0.7051    4.7646    1.3479    9.3722   -1.4222    7.9275
% 7.5777    1.7549    1.1836]
disp(sum(randfac))
% E = 38.7019
John Alexiou
  • 28,472
  • 11
  • 77
  • 133
  • Wow, yeah, that might be ideal. Better than what I had for sure. Let me think about it for a bit, and then I'll probably accept it. Thanks a lot. – user129412 Dec 06 '14 at 17:29
  • It's worth noting that this doesn't sample from a uniform distribution. It's closer to a normal distribution. – Veedrac Dec 06 '14 at 19:20
1

Just divide the vector with the sum and multiply with the target E.

randfac = (maxrand-minrand).*rand(1,N) + minrand;
randfac = E*randfac/sum(randfac);

as long as the operator is linear, the result is going to retain it's randomness. Below is some sample code:

minrand = 0;
maxrand = 1;
N = 1000; %size
v = (maxrand-minrand).*rand(1,N) + minrand;
E = 100; %Target sum
A = sum(v);
randfac = (E/A)*v;

disp(sum(randfac))
% 100.0000
John Alexiou
  • 28,472
  • 11
  • 77
  • 133
  • I don't see this. We have the vector V, for which the sum is equal to E. So we define v = V/E. We then multiply this by the random numbers multiplied by E. That just cancels the division by E, giving us V*random vector. There's no guarantee that this will sum to E, right? – user129412 Dec 06 '14 at 16:22
  • If the sum of a vector `v` is `A` then the `sum(2*v)=2*A` because sum is a linear operator. So if `sum(v)=A` then `sum(v/A)=1` and thus `sum(E*(v/A))=E` – John Alexiou Dec 06 '14 at 16:24
  • Ahh, I'm stupid. I missed the /sum(randfac), so I didn't realize it had to sum to E as well. Hm, I guess this could work! – user129412 Dec 06 '14 at 16:30
  • Yes, it does work. Consider the sum `x=a+b+c` then what is `2a+2b+2c`? It is `2x`. Do you see it now? – John Alexiou Dec 06 '14 at 16:31
  • Hm, well I guess I wasn't clear in what I wrote in my original post. the vector V is given; I don't generate it. It sums to E. After that, I want to do the randomization stuff, making sure it still sums to E. I'm afraid you slightly misunderstood? – user129412 Dec 06 '14 at 16:35
  • The way I understood it you simply mean e = e/E; minrand = 0; maxrand = 2; randfac = (maxrand-minrand).*rand(1,N) + minrand; randfac = E*randfac/sum(randfac); e = e.*randfac; – user129412 Dec 06 '14 at 16:36
  • Please update the question if you have not explained it well. I have answered per the original question. – John Alexiou Dec 06 '14 at 17:09
1

First of all with random numbers in the interval of [a b] you can't guarantee that you will have the same summation (same E). For example if [a b]=[1 2] of course the E will increase.

Here is an idea, I don't know how random is this!

For even N I randomize V then divide it in two rows and multiply one of them with random numbers in [a b] but the second column will be multiplied to a vector to hold the summation fixed.

N = 10;
V = randi(100,[1 N]);
E = sum(V);
idx = randperm(N);
Vr = V(idx);
[~,ridx] = sort(idx);
Vr = reshape(Vr,[2 N/2]);
a = 1;
b = 3;
r1 = (b - a).*rand(1,N/2) + a;
r2 = (sum(Vr) - r1.*Vr(1,:))./Vr(2,:);
r = reshape([r1;r2],1,[]);
r = r(ridx);
Enew = sum(V.*r);

The example results are,

V = [12      82      25      51      81      51      31      87      6       74];
r = [2.8018  0.7363  1.9281  0.5451  1.9387 -0.4909  1.3076  0.8904  2.9236  0.8440];

with E = 500 as well as Enew.

I'm simply assigning one random number to a pair (It can be considered as half random!).

Rashid
  • 4,326
  • 2
  • 29
  • 54
  • You are right, I would have to choose [a,b] such that (b-a)/2 is smaller or equal to 1. – user129412 Dec 06 '14 at 17:12
  • I like this approach, however for my purposes it might not work too well because it doesn't treat the first N/2 elements the same way as it does the second N/2, and this is a bit of an issue. However, I didn't write that in my problem, and this is definitely a solution to the way I have stated it! – user129412 Dec 06 '14 at 17:14
  • @user129412 Yes, But if you first randomize `V` then apply it, it will be random enough! let me update the post. – Rashid Dec 06 '14 at 17:15
  • @user129412, I edited the post. I think this is the most random that my solution can get! I hope it helps, Good luck. – Rashid Dec 06 '14 at 17:29
0

Okay, I have found a way to somewhat do this, but it is not elegant and there are probably better solutions. Starting with an initial vector e, for which sum(e) = E, I can randomize its values and end up with an e for which sum(e) is in the range [(1-threshold)E,(1+thresholdE)]. It is computationally expensive, and not pretty.

The idea is to first multiply e by random numbers in a certain range. Then, I will check what the sum is. If it is too big, I will decrease the value of the random numbers smaller than half of the range until the sum is no longer too big. If it is too small, I do the converse, and iterate until the sum is within the desired range.

e = somepredefinedvector
minrand = 0;
maxrand = 2;
randfac = (maxrand-minrand).*rand(1,N) + minrand;
e = randfac.*e;
threshold = 0.001;

while sum(e) < (1-threshold)*E || sum(e) > (1+threshold)*E
    if sum(e) > (1+threshold)*E
        for j = 1:N
            if randfac(j) > (maxrand-minrand)/2
                e(j) = e(j)/randfac(j);
                randfac(j) = ((maxrand-minrand)/2-minrand).*rand(1,1) + minrand;
                e(j) = randfac(j)*e(j);
            end
            if sum(e) > (1-threshold)*E && sum(e) < (1+threshold)*E
                break
            end
        end
    elseif sum(e) < (1-threshold)*E
        for j = 1:N
            if randfac(j) < (maxrand-minrand)/2
                e(j) = e(j)/randfac(j);
                randfac(j) = (maxrand-(maxrand-minrand)/2).*rand(1,1) + (maxrand-minrand)/2;
                e(j) = randfac(j)*e(j);
            end
            if sum(e) > (1-threshold)*E && sum(e) < (1+threshold)*E
                break
            end
        end
    end
end
user129412
  • 765
  • 1
  • 6
  • 19