3
\$\begingroup\$

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?

\$\endgroup\$
  • \$\begingroup\$ There's no need to state what kind of improvements you're looking for in the title. The moment you post your code on Code Review it's fair game for the reviewers to critique any and all parts of the code. \$\endgroup\$ – Mast Nov 16 '15 at 15:47
  • \$\begingroup\$ Ok! For specific questions I will post specific lines code then! Thanks! \$\endgroup\$ – John Ludos Nov 17 '15 at 15:04
1
\$\begingroup\$

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.

\$\endgroup\$
  • \$\begingroup\$ Thanks, good idea to set out the construction of the cumulative weight. Any idea about how to vectorise the for cycle? This may not change significantly the computational time... but I'm curious how this can be done technically! \$\endgroup\$ – John Ludos Nov 17 '15 at 15:06
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.