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;    endenddisp(beans)Any idea how to vectorise the two for-cycle in the function?
Weighted dice with n faces
matlab;dice;vectorization
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));endwhich 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   endend% assign the interval found to the value of the faceface_index = j;end%%%--------------main scriptcumulative_weight = zeros(1, num_weight + 1);for j=1:num_weight   cumulative_weight(j+1) = sum(weight(1:j));endval = zeros(1,1000);for i=1:1000    val(i) = new_weighted_dice(cumulative_weight)endbeans = histc(val, cumulative_weight);Short, clean, fast.