Probability of binary sequence with Monte Carlo

6 vues (au cours des 30 derniers jours)
pepito che
pepito che le 20 Nov 2020
Commenté : pepito che le 23 Nov 2020
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

Réponse acceptée

James Tursa
James Tursa le 20 Nov 2020
Modifié(e) : James Tursa le 20 Nov 2020
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
  3 commentaires
James Tursa
James Tursa le 22 Nov 2020
Modifié(e) : James Tursa le 22 Nov 2020
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
pepito che le 23 Nov 2020
Yes I get it. Thanks for your help :)

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur PHY Components dans Help Center et File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by