Probability of binary sequence with Monte Carlo

Hello everyone !

Context:

I'm trying to compute the probability that a binary sequence "stop" under a certain constraint: if the sequence has a majority of "1" we stop.

ex: Let's say n is the length of the sequence; At n=2 we have : {(00),(01),(10),(11)} possibilities and so at n=3 {(000),(001),(010),(011),(100),(101)} possibilities (note: we didn't took (1,1) car we already stop at n=2) and so on... the probability to stop at n=2 is P(2)=1/4 and at n=3, P(3)=2/6

Question:

My goal is to know what is the probability that at "step n" we stop (ie, how many sequence at n-step "stop" over how many possibilities, taking in account tthat some possibility already stop before)?

Tentative:

The probleme is the more n will grow up the bigger possibilities we have and im not sure the simulation can support for n tend to big number.

So I heard about "Monte carlo simulation"; my idea was to reapeat a certain time the experience. Like I generate a random binary sequence of length-n and i look the majority (ie if sum(sequence)>=n/2 , it's important that this appear in code) and try to generate statistic. I don't know how much "monte carlo simulation" could help me in my problem.

I did this:

rng('shuffle')

nb_mort=0;  %number of sequence who stop
n=0;       %length of sequence
maj=0;     %majority
M=1000;  %nb of run


for n_run=1:M;
    
    while (n <= 100) | (maj <= n/2 )
        
        
        count=zeros(n);
        n=n+1
        
        seq = randi([0 1],1,n)
        
        maj=sum(seq)
                       
    end
    
    clear n;    %I want he start again the loop but renitialize his length
    n=0;
    
       nb_mort=nb_mort+1;
                       
    %count(n)=count(n)+1;
    
    %proba_n(n)=count(n)/M;
    
    %nb_mort_mean(n_run)=nb_mort
    

end

Don't take in account, the last line in comment Proba, because the proba is wrong, i didn't found...

I don't know why also at the end of my "while" he doesnt renialize the n??

Im not super familiar with Monte Carlo and im beginner on Matlab

Thank you for your help

Answer 1

[James Tursa]

First point is that probability of stopping at n=3 is not 2/6, it is (3/4)*(2/6) because you have to include the probablility that you didn't stop at n=2. Regardless ...

Seems like you could simply generate a long string of m digits and then use cumsum on that compared to 1:m to see where you stopped for that string. I.e., first index where the cumsum element is greater than the (1:m)/2 values is the stopping point. If you didn't stop in m digits put all of those results into a single "greater than m" category. Wrap that in a large Monte-Carlo loop and then divide your counts by the number of Monte-Carlo runs to get the individual probabilities.

And it looks like you have a rule to never stop at n=1. True?

E.g., initialize some values

m = 100; % some value for maximum stopping point we will look for
count = zeros(1,m+1);
m2 = (1:m)/2; % the comparison vector

And then each iteration of the Monte-Carlo loop would be something like

seq = randi([0 1],1,m);
cseq = cumsum(seq);
cseq(1) = 0; % force that we won't stop at n=1
x = find(cseq > m2, 1); % find the first stopping point
if isempty(x)
    count(m+1) = count(m+1) + 1;  % we didn't stop, so lump this into special category
else
    count(x) = count(x) + 1; % increment the 1st place we stopped
end

[pepito che]

Thank you for your help !

I didn't know the command "cumsum" which indeed is very usefull for my problem.

Indeed i cannot stop at n=1 because i always suppose that i start minimum with a sequence of 2-bit.

Just a point where Im not completlely clear is the:

count(m+1) = count(m+1) + 1;  % we didn't stop, so lump this into special category

In here you put at the end of the vector "count", the number of time when we didn't stop; right ?

In theroy i could remove this line no ? If im only interesting about the case when "we stop".

[James Tursa]

No matter what value of m you select, there will always be a probability that you didn't stop in m bits. In general you should expect that a certain percentage of your iterations will be this case. So you need to account for those cases. So yes, the count(m+1) spot counts all the times that we didn't stop. If you don't care about that then of course you can ignore that part of the result and not even calculate it, but it will be a non-zero value so you need logic to detect the "x is empty" case in some way. E.g., if you didn't even care to calculate this probability then you could do this:

if ~isempty(x)
    count(x) = count(x) + 1; % increment the 1st place we stopped
end

That being said, I would advise you to not do this. These cases actually happen, so why not calculate the probability since you already have all the logic for it?

[pepito che]

Yes I get it. Thanks for your help :)