Weighted dice with n faces

Question

Here you find some code that returns the value of a multivariate Bernoulli density function - for a dice with n faces, where each face has its weight (weight must be positive reals, not necessarily sum = 1) it provides the result of the tossing.

function face = weighted_dice(faces, weight)
    % face = dice(faces, weight) 
    % input: row vector. Faces and weight of the same dimension!
    % output: face is the discrete face after tossing a virtual dice with given
    % faces and weight.
    if size(faces, 2) ~= size(weight, 2) || size(weight, 2)<= 1
        error('Input of function dice not well defined. See help')
    end

    num_weight = size(weight, 2);
    % weight = [w1,w2,w3, .. ,w_n]
    % cumulative_weight = [0, w1, w1+w2, w1+w2+w3, ... , sum(weight)]

    cumulative_weight = zeros(1, num_weight + 1);
    for j=1:num_weight
       cumulative_weight(j+1) = sum(weight(1:j));
    end

    % pick a random value unif. distrib. betw. 0 and the sum of the weights.
    rand_w = unifrnd(0, cumulative_weight(end));
    % check in which interval of the cumulative_weight is this element 
    for j=1:num_weight
       if rand_w >= cumulative_weight(j) && rand_w < cumulative_weight(j+1)
           break
       end
    end

    % assign the interval found to the value of the face
    face = faces(j);

E.g.

disp(weighted_dice(['a','b','c','d','e','f'], [0.1, 0.1, 0.5, 0.1, 0.1, 0.1]))

beans = zeros(1,6);
for i=1:1000
    val = weighted_dice(['a','b','c','d','e','f'], [0.1, 0.1, 0.5, 0.1, 0.1, 0.1]);
    if val == 'a'
        beans(1,1) = beans(1,1) + 1;
    elseif val == 'b'
        beans(1,2) = beans(1,2) + 1;
    elseif val == 'c'
        beans(1,3) = beans(1,3) + 1;
    elseif val == 'd'
        beans(1,4) = beans(1,4) + 1;
    elseif val == 'e'
        beans(1,5) = beans(1,5) + 1;
    elseif val == 'f'
        beans(1,6) = beans(1,6) + 1;
    end
end

disp(beans)

Any idea how to vectorise the two for-cycle in the function?

Answer 1

If you're trying to speed up the code, I would do a couple of things. You have some redundant processing, such as this:

cumulative_weight = zeros(1, num_weight + 1);
for j=1:num_weight
   cumulative_weight(j+1) = sum(weight(1:j));
end

which is performed 1000 times inside weighted_dice, but you really only need to do it once. I would also get rid of the face labels altogether, because you don't care what the 6 faces are called, you just care about the distribution of results. Finally, you can use the MATLAB builtin histc to generate the results in a vectorized fashion.

function face_index = new_weighted_dice(cumulative_weight)
    % cumulative_weight = cumulative distribution function to draw against
    % output: face_index is the discrete face after tossing a virtual dice with given
    % faces and weight.


% pick a random value unif. distrib. betw. 0 and the sum of the weights.
rand_w = unifrnd(0, cumulative_weight(end));
% check in which interval of the cumulative_weight is this element 
for j=1:num_weight
   if rand_w >= cumulative_weight(j) && rand_w < cumulative_weight(j+1)
       break
   end
end

% assign the interval found to the value of the face
face_index = j;
end

%%%--------------main script

cumulative_weight = zeros(1, num_weight + 1);
for j=1:num_weight
   cumulative_weight(j+1) = sum(weight(1:j));
end
val = zeros(1,1000);
for i=1:1000
    val(i) = new_weighted_dice(cumulative_weight)
end

beans = histc(val, cumulative_weight);

Short, clean, fast.

Answer 2

If you are willing to give up having an actual weighted_dice() function, then a very simple way to generate the data you want is

beans = histc(cumulative_weight, rand(1000,1))

Note: this generates a 7x1 array, and you only want the first 6 entries. Also, I don't have the stats toolbox on this machine, so I used rand, but since you're working on the interval [0,1] they are equivalent.