Search code examples
matlabbinomial-cdf

Code returns minimal m s.t P(X>m)=1/n^2 for binomial experiment


Assume we do an experiment where we throw n balls into n jars. X is an independent variable describes number of balls int the first jar. Build a function returns the smallest integer fulfills P(X>m)<1/n^2.

The distribution is binomial so I wrote the following matlab function:

function m = intpq3(n)
flag=0;
par=1/n^2;
m=0;
P=0;
%Since X is non-neative integer 
if(n==1)
    m=-1*Inf;
else
    while(flag==0 && m<=n)
        m=m+1;
         P=P+nchoosek(n,m)*(1/n)^m*(1-1/n)^(n-m);   

        if(1-P<=par)
            flag=1;
        end

    end
end
disp(m)
end

But for every 'n' I give it, it returns either error or n-1. What am I doing wrong?


Solution

  • The following version of your program seems to do what you want. The problem with your version as far as I can tell is that you did not include m=0 into your sum of P, thus 1-P was consistently too large. It's always a good idea to ask the program to spit out numbers and compare to paper-and-pencil calculations.

    function [m,P] = intpq3(n,varargin);    
    
    if nargin==2
       plt=varargin{1}; 
    else
       plt = 0;
    end
    
    flag=0;
    par=1/n^2;
    m=0;
    P=0;
    
    disp(['1/n^2 =  ' num2str(par)]);
    
    %Since X is non-neative integer 
    if(n==1)
        m=-1*Inf;
    else
        while m<n & flag==0
           P=P+nchoosek(n,m)*(1/n)^m*(1-1/n)^(n-m);   
           disp([' n= ' num2str(n)  ' m= ' num2str(m)  ' P(X>' num2str(m) ') ' num2str(1-P)])
            if(1-P<=par)
                flag=1;
            end
           m=m+1;
        end
    end
    disp(['mselect =  ' num2str(m-1)]);
    end
    

    The following also helped to diagnose the problem. It shows values of P(X>m) (colored dots) for different m at select values of n. Overlayed as a dashed line is 1/n^2.

    enter image description here