matlabperformancestatisticsprobabilitybinomial-cdf

Calculating the probability of success in k (or less) Bernoulli trials out of n using matlab


I am trying to calculate the probability of success in 70 (or less) Bernoulli trials out of 100. I wrote it with Matlab. But, I get the probability to be 1 (it can't be 1 since its not success in all 100 trials).

Is my function OK?

syms k
f = nchoosek(100,k)*0.5^k*0.5^(100-k);
F = double(symsum(nchoosek(100,k)*0.5^k*0.5^(100-k),k,0,70));

If it is, how can I get a more accurate resault in Matlab?

Thanks

edit: I have have a binary vector that represents success/failure in n trials (like tossing a coin 100 times). And I need the error of my sample (the way statistics does it.. but I don't know statistics). So I thought that maybe i will try to calculate "how far am I from being correct in all trials" which should be 1-F in my code. But then 70 successes out of 100 gives me error = 0 which is obviously not true..

edit2: In the example I gave here I need the probability that there are 70 successes in 100 trials.

enter image description here


Solution

  • You do have everything you need to answer this question.

    In the formula you have posted, you sum the probabilities from 0 to 70, that is, it will calculate the probability to have 0 or 1 or 2 .. or 70 successes, which means 70 or less successes.

    Without the sum, you get the probability to have exactly k successes. The probability to get exactly 70 successes is:

    k = 70;
    f = nchoosek(100,k)*0.5^k*0.5^(100-k)
    Warning: Result may not be exact. Coefficient is greater than 9.007199e+15 and is only
    accurate to 15 digits 
    > In nchoosek (line 92) 
    
    f =
    
       2.3171e-05
    

    You receive a warning that the computation of nchoosek(100,70) is not exact (see below for a better way).

    To compute the probability to get 70 or less successes, sum over the probabilities to get 0 or 1 or .. 70 successes:

    >> f = 0;
    >> for k=0:70;
    f = f + nchoosek(100,k)*.5^k*.5^(100-k);
    end
    

    You will receive a lot of warnings, but you can look at f:

    >> f
    
    f =
    
        1.0000
    

    As you see, if rounded to four digits, the probability is 1. We know, however, that it must be slighly less than oneĀ“. If we ask Matlab to show more digits:

    >> format long
    

    we see that it is not exactly 1:

    >> f
    
    f =
    
       0.999983919992352
    

    If you compute 1-f, you will see that the result is not 0 (I switch back to showing less digits):

    >> format short
    >> 1-f
    
    ans =
    
       1.6080e-05
    

    To get rid of the warnings and to simplify the code for computing the probabilities, Matlab provides several functions to deal with binomial distributions. For the probability to get exactly 70 successes, use

    >> binopdf(70,100,.5)
    
    ans =
    
       2.3171e-05
    

    and to get 70 or less successes:

    >> format long
    >> binocdf(70,100,.5)
    
    ans =
    
       0.999983919992352