Small Probabilities appearing less frequent using rand

I am trying to do some Monte Carlo Type simulations. In this, I have few choices which needs to be selected using weighted probability. To do this, I have calculated a cumulative probability (Prob_Rv) and I access the choice using something like

choice = find(Prob_Rv > rand, 1, 'first');

For a example scenario, Lets take Prob_Rv as

Prob_Rv = [0.0001;0.0001;0.9993;1.0000];

From this example, the Probability of choice 2 is zero, however rest of the 3 choices can be with certain probability. On doing something like this, I find that some of the choice occur more frequently than they should and some less frequently than they should. I tried to test it using a code like this:

    Prob4 = 7.4831e-4;
    totNum = 2000000;
    nn = zeros(totNum,1);
    storeVal = zeros(100,4);
    for i = 1:100
        rng shuffle;
        nn = rand(totNum,1);
        Prob1_1 = (sum(nn < Prob_Rv(1)))/totNum;
        Prob2_1 = (sum(nn < Prob_Rv(2))/totNum) - Prob1_1;
        Prob3_1 = (sum(nn < Prob_Rv(3))/totNum) - Prob2_1;
        Prob4_1 = 1 - Prob1_1 - Prob2_1 - Prob3_1;
        storeVal(i,:) = [Prob1_1 Prob2_1 Prob3_1 Prob4_1];
    end
    plot(storeVal(:,4));
    hold on;
    plot(1:100,Prob4*ones(100,1),'*');

In this, I create 200,000 random numbers and check the probability of a choice. In the image below, I have plotted the probability obtained for choice 4 and compared it with the probability it should be. The probability I get are nearly always lower than it should be.

Is it something I am doing wrong? Is there a way to resolve this?

Edit: #1 I have tried removing the 'rng shuffle' in every loop, and the problem remains the same.

[Amit]

This was my mistake. Prob3_1 should have been

Prob3_1 = (sum(nn < Prob_Rv(3))/totNum) - Prob2_1 - Prob1_1;

I don't know how to delete the question otherwise I would have.

Answer 1

[the cyclist]

Your code gives an error because Prob4 is undefined. I changed that variable to Prob4_1, and the random values are well centered on the expected value.

[Amit]

The code does not gives error. The Prob4 is actual probability (and defined earlier. I only posted a portion of the code)and plotted for comparison (symbol '*').

From the figure it can be seen that '*' (actual probability) is approximately 1.2 times higher than given by (on average) matlab.

Answer 2

[Roger Stafford]

Your computation is faulty. Where you write

    Prob2_1 = (sum(nn < Prob_Rv(2)) - Prob1_1)/totNum;
    Prob3_1 = (sum(nn < Prob_Rv(3)) - Prob2_1)/totNum;

it should be:

    Prob2_1 = sum(nn < Prob_Rv(2))/totNum - Prob1_1;
    Prob3_1 = sum(nn < Prob_Rv(3))/totNum - Prob2_1;

You already divided by 'totNum' to get 'Prob2_1'. You shouldn't be doing a second time in 'Prob3_1'. The same applies for division by 'totNum' in 'Prob1_1' and 'Prob2_1', though in your case it makes no difference.

That may account for the discrepancy you describe.

[Amit]

I am sorry, I think I posted the code from second try (wrong code). I have updated the code as well as the plot obtained.

The issue still exists.

[Roger Stafford]

Please post the exact code and values for which the issue still exists. Please give the values using "format long" since you are dealing in probabilities of such a small magnitude.