[Contents] [TitleIndex] [WordIndex

Source code of Binomial macros

Note -- latex needs to be added , and some more examples also need to be added .

Introduction

This is a report by Prateek Papriwal for the GSOC 2012 on the project "Distribution functions". This page gathers proposals for the Binomial distribution.

A word of caution

These function are just temporary work and may contain bugs, which may have been fixed in the distfun project at http://forge.scilab.org/index.php/p/distfun/source/.

distfun_binocdf.sci

// Copyright (C) 2012 - Prateek Papriwal
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//

function p = distfun_binocdf(varargin)
    // Binomial CDF
    //
    // Calling Sequence
    //   p = distfun_binocdf(X,N,Pr)
    //   p = distfun_binocdf(X,N,Pr,lowertail)
    //
    //Parameters
    //   X : a 1x1 or nxm matrix of doubles, the number of Bernoulli trials after in which success occurs . X belongs to the set {0,1,2,3,.....}
    //   N : a 1x1 or nxm matrix of doubles , the total number of binomial trials . N belongs to the set {1,2,3,4,.......}
    //   Pr : a 1x1 or nxm matrix of doubles,  the probability of success in a Bernoulli trial
    //   lowertail : a 1x1 matrix of booleans, the tail (default lowertail=%t). If lowertail is true (the default), then considers P(X<x) otherwise P(X>x).
    //   p : a nxm matrix of doubles, the probability.
    //
    // Description
    //   Computes the cumulative distribution function of 
    //   the Binomial distribution function.
    //   
    //   Any scalar input argument is expanded to a matrix of doubles 
    //   of the same size as the other input arguments.
    //
    //Examples
    // // Check with X scalar, N scalar, Pr scalar
    // computed = distfun_binocdf(100,162,0.5);
    // expected = 0.9989567;
    //
    // // Check with expanded X
    // computed = distfun_binocdf([5 15],100,0.05);
    // expected = [0.6159991 0.9999629];
    //
    // // Check with expanded N
    // computed = distfun_binocdf(5,[50 100],0.05);
    // expected = [0.9622238 0.6159991];
    //
    // // Check with expanded Pr
    // computed = distfun_binocdf(5,50,[0.05 0.1]);
    // expected = [0.9622238 0.6161230];
    //
    // // Check with two arguments expanded 
    // computed = distfun_binocdf([5 10],[50 100],0.05);
    // expected = [0.9622238 0.9885276];
    //
    // // Check with all the arguments expanded
    // computed = distfun_binocdf([5 10],[50 100],[0.05 0.1]);
    // expected = [0.9622238 0.5831555];
    //
    // //Plot the function
    //
    // scf();
    // Pr1 = 0.5;
    // Pr2 = 0.7;
    // Pr3 = 0.5;
    // for i = 0:40
    //     X = linspace(i,i+1,20);
    //     y = distfun_binocdf(i,20,Pr1);
    //     plot(X,y,'b');
    //     y = distfun_binocdf(i,20,Pr2);
    //     plot(X,y,'g');
    //     y = distfun_binocdf(i,40,Pr3);
    //     plot(X,y,'r');
    // end
    // xtitle("Geometric CDF","x","P(X<x)");
    //
    // //check upper tail
    //p = distfun_binocdf(3,10,0.1);
    //lt_expected = 0.9872048;
    //
    //q = distfun_binocdf(3,10,0.1,%f);
    //ut_expected = 0.0127952;
    //
    // Bibliography
    // http://en.wikipedia.org/wiki/Binomial_distribution
    //
    // Authors
    // Copyright (C) 2012 - Prateek Papriwal
    
    [lhs,rhs] = argn()
    apifun_checkrhs("distfun_binocdf",rhs,3:4)
    apifun_checklhs("distfun_binocdf",lhs,0:1)
    
    X = varargin(1)
    N = varargin(2)
    Pr = varargin(3)
    lowertail = apifun_argindefault(varargin,4,%t)
    //
    // Check type
    //
    apifun_checktype("distfun_binocdf",X,"X",1,"constant")
    apifun_checktype("distfun_binocdf",N,"N",2,"constant")
    apifun_checktype("distfun_binocdf",Pr,"Pr",3,"constant")
    apifun_checktype("distfun_geocdf",lowertail,"lowertail",4,"boolean")
    
    apifun_checkscalar("distfun_geocdf",lowertail,"lowertail",4)
    //
    //Check content
    // 
    apifun_checkgreq("distfun_binocdf",X,"X",1,0)
    apifun_checkgreq("distfun_binocdf",N,"N",2,1)
    apifun_checkrange("distfun_binocdf",Pr,"Pr",3,0,1)
    
    [X,N,Pr] = apifun_expandvar(X,N,Pr)
    nx = size(X,"c")
    p = ones(1,nx) - 1
    for j = 1:nx 
        for i = 0:X(j) 
            p(1,j) = p(1,j) + distfun_binopdf(i,N(j),Pr(j)) 
        end
    end
    if (X == []) then
        p=[]
        return
    elseif (N == []) then
        p = []
    elseif (Pr == []) then
        p=[]
        return
    elseif (lowertail == %f) then
        path = distfun_getpath()
        internallib = lib(fullfile(path,"macros","internals"))
        p = distfun_p2q(p)
    else
    
    end
endfunction

distfun_binoinv.sci

// Copyright (C) 2012 - Prateek Papriwal
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
//

function X = distfun_binoinv(varargin)
    // Binomial Inverse CDF
    //
    // Calling Sequence
    //   X = distfun_binocdf(p,N,Pr)
    //   X = distfun_binocdf(p,N,Pr,lowertail)
    //
    //Parameters
    //   p : a nxm matrix of doubles, the probability.
    //   N : a 1x1 or nxm matrix of doubles , the total number of binomial trials . N belongs to the set {1,2,3,4,.......}
    //   Pr : a 1x1 or nxm matrix of doubles,  the probability of success in a Bernoulli trial
    //   lowertail : a 1x1 matrix of booleans, the tail (default lowertail=%t). If lowertail is true (the default), then considers P(X<x) otherwise P(X>x).
    //   X : a 1x1 or nxm matrix of doubles, the number of Bernoulli trials after in which success occurs . X belongs to the set {0,1,2,3,.....}
    //
    // Description
    //   Computes the Inverse cumulative distribution function of 
    //   the Binomial distribution function.
    //   
    //   Any scalar input argument is expanded to a matrix of doubles 
    //   of the same size as the other input arguments.
    //
    //Examples
    //
    // // Check with all arguments scalar
    // X = distfun_binoinv(0.21,162,0.5);
    // expected = 76;
    //
    // // check with p expanded
    // X = distfun_binoinv([0.05 0.95],162,0.5);
    // expected = [71 91];
    //
    // // Check with N expanded
    // X = distfun_binoinv(0.05,[100 162],0.5);
    // expected = [42 71];
    //
    // // Check with expanded Pr
    // X = distfun_binoinv(0.05,162,[0.2 0.5]);
    // expected = [24 71];
    //
    // //Check with all arguments expanded
    // X = distfun_binoinv([0.05 0.95],[100 162],[0.2 0.5]);
    // expected = [14 91];
    //
    // Bibliography
    // http://en.wikipedia.org/wiki/Binomial_distribution
    //
    // Authors
    // Copyright (C) 2012 - Prateek Papriwal

    [lhs,rhs] = argn()
    apifun_checkrhs("distfun_binoinv",rhs,3:4)
    apifun_checklhs("distfun_binoinv",lhs,0:1)
    
    p = varargin(1)
    N = varargin(2)
    Pr = varargin(3)
    lowertail = apifun_argindefault(varargin,4,%t)
    //
    //Check type
    //
    apifun_checktype("distfun_binoinv",p,"p",1,"constant")
    apifun_checktype("distfun_binoinv",N,"N",2,"constant")
    apifun_checktype("distfun_binoinv",Pr,"Pr",3,"constant")
    apifun_checktype("distfun_binoinv",lowertail,"lowertail",4,"boolean")


    apifun_checkscalar("distfun_binoinv",lowertail,"lowertail",4)
    //
    //Check Content
    //
    apifun_checkrange("distfun_binoinv",p,"p",1,0,1)
    apifun_checkgreq("distfun_binoinv",N,"N",2,1)
    apifun_checkrange("distfun_binoinv",Pr,"Pr",3,0,1)
    
    [p,N,Pr] = apifun_expandvar(p,N,Pr)
    
    path = distfun_getpath()
    internallib = lib(fullfile(path,"macros","internals"))
    q = distfun_p2q(p)
    
    if(p == []) then
        X = []
        return
    elseif(N == []) then
        X = []
        return
    elseif(Pr == []) then
        X = []
        return
    elseif(lowertail) then
        X = ceil(distfun_invcdfbin(N,Pr,1-Pr,p,1-p))
    else
        X = floor(distfun_invcdfbin(N,Pr,1-Pr,1-p,p))
    end
endfunction

distfun_binopdf.sci

// Copyright (C) 2012 - Prateek Papriwal
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function y = distfun_binopdf(varargin)
    // Binomial PDF
    //
    //Calling Sequence
    //y = distfun_binopdf(X,N,Pr) 
    //
    //Parameters
    //   X : a 1x1 or nxm matrix of doubles, the number of Bernoulli trials after in which success occurs . X belongs to the set {0,1,2,3,.....}
    //   N : a 1x1 or nxm matrix of doubles , the total number of binomial trials . N belongs to the set {1,2,3,4,.......}
    //   Pr : a 1x1 or nxm matrix of doubles,  the probability of success in a Bernoulli trial 
    //   y : a nxm matrix of doubles, the Probability density.
    //
    // Description
    //   Computes the probability distribution function of 
    //   the Binomial distribution function.
    //   
    //   Any scalar input argument is expanded to a matrix of doubles 
    //   of the same size as the other input arguments.
    
    
    //Examples
    // // Check with X scalar, N scalar, Pr scalar
    //y = distfun_binopdf(0,200,0.02);
    //expected = 0.0175879;
    
    // // Check with expanded X
    //computed = distfun_binopdf([5 15],100,0.05);
    //expected = [0.1800178 0.0000988];
    
    // // Check with expanded N
    //computed = distfun_binopdf(5,[50 100],0.05);
    //expected = [0.0658406 0.1800178];
    
    // // Check with two arguments expanded 
    //computed = distfun_binopdf([5 10],[50 100],0.05);
    //expected = [0.0658406 0.0167159];
    
    // // Check with all the arguments expanded
    //computed = distfun_binopdf([5 10],[50 100],[0.05 0.1]);
    //expected = [0.0658406 0.1318653];
    
    // // Plot the function
    // // Note that geometric distribution is discrete in nature 
    // // i.e it gives output only at discrete values of X = {0,1,2,3,.....}
    
    // scf();
    // N1 = 20;
    // X = 0:N1
    // y1 = distfun_binopdf(X,N1,0.5);
    // plot(X,y1,"bo-");
    // N2 = 20;
    // X = 0:N2
    // y2 = distfun_binopdf(X,N2,0.7);
    // plot(X,y2,"go-");
    // N3 = 40;
    // X = 0:N3;
    // y3 = distfun_binopdf(X,N3,0.5)
    // plot(X,y3,"ro-");
    
    // Bibliography
    // http://en.wikipedia.org/wiki/Binomial_distribution
    //
    // Authors
    // Copyright (C) 2012 - Prateek Papriwal
    
    [lhs,rhs] = argn()
    apifun_checkrhs("distfun_binopdf",rhs,3)
    apifun_checklhs("distfun_binopdf",lhs,0:1)
    
    X = varargin(1)
    N = varargin(2)
    Pr = varargin(3)
    
    //
    // Check type
    //
    apifun_checktype("distfun_binopdf",X,"X",1,"constant")
    apifun_checktype("distfun_binopdf",N,"N",2,"constant")
    apifun_checktype("distfun_binopdf",Pr,"P",3,"constant")
    //
    // Check content
    //  
    apifun_checkgreq("distfun_binopdf",X,"X",1,0)
    apifun_checkgreq("distfun_binopdf",N,"N",2,1)
    apifun_checkrange("distfun_binopdf",Pr,"P",3,0,1)
    
    
    [X,N,Pr] = apifun_expandvar(X,N,Pr)
    
    y = %e .^ (gammaln(N+1)-gammaln(X+1)-gammaln(N-X+1)) .* (Pr .^ X) .* ((1-Pr) .^(N-X))
endfunction

distfun_binornd.sci

// Copyright (C) 2012 - Prateek Papriwal
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function R = distfun_binornd(varargin)
    // Binomial random numbers
    //
    // Calling Sequence
    //   R = distfun_binornd(N,Pr)
    //   R = distfun_binornd(N,Pr,v)
    //   R = distfun_binornd(N,Pr,m,n)
    //
    // Parameters
    //   N : a 1x1 or nxm matrix of doubles , the total number of binomial trials . N belongs to the set {1,2,3,4,.......}
    //   Pr : a 1x1 or nxm matrix of doubles, the probability of getting success in a Bernoulli trial
    //   v : a 1x2 or 2x1 matrix of doubles, the size of R
    //   v(1) : the number of rows of R
    //   v(2) : the number of columns of R
    //   m : a 1x1 matrix of floating point integers, the number of rows of R
    //   n : a 1x1 matrix of floating point integers, the number of columns of R
    //   R: a matrix of doubles, the random numbers.
    //
    //
    // Description
    //   Generates random variables from the Binomial distribution function.
    //
    //   Any scalar input argument is expanded to a matrix of doubles 
    //   of the same size as the other input arguments.
    //
    //   As a side effect, it modifies the internal seed of the grand function.
    //
    // Note - The output argument R belongs to the set {0,1,2,3,...}.
    //
    // Examples
    //
    // // set the initial seed for tests
    // distfun_seedset(1);
    
    // // Check with expanded N
    // N = [10 100 1000 10000];
    // Pr = 0.1;
    // computed = distfun_binornd(N, Pr);
    // expected = [1 19 98 964];
    //
    // // set the initial seed
    // distfun_seedset(1);
    // // Check with expanded Pr
    // N =100;
    // Pr = [0.1 0.2 0.3 0.4];
    // computed = distfun_binornd(N, Pr);
    // expected = [9 32 29 38];
    //
    // //Check R = distfun_binornd(Pr,v)
    // computed = distfun_binornd(100,0.2,[4 5]);
    // assert_checkequal(size(computed),[4 5]);
    //
    // //Check mean and variance
    // N = 1000;
    // Pr = 0.3;
    // n = 5000;
    // computed = distfun_binornd(N,Pr,[1 n]);
    // c = mean(computed(1:n));
    // d = st_deviation(computed(1:n) );
    // [M,V] = distfun_geostat (Pr);
    //
    // Bibliography
    // http://en.wikipedia.org/wiki/Binomial_distribution
    //
    // Authors
    // Copyright (C) 2012 - Prateek Papriwal
    //
    path = distfun_getpath()
    internallib = lib(fullfile(path,"macros","internals"))
    
    [lhs,rhs] = argn()
    apifun_checkrhs("distfun_binornd",rhs,2:4)
    apifun_checklhs("distfun_binornd",lhs,0:1)
    
    N = varargin(1)
    Pr = varargin(2)
    //
    // Check type
    //
    apifun_checktype("distfun_binopdf",N,"N",1,"constant")
    apifun_checktype("distfun_binornd",Pr,"Pr",2,"constant")
    
    if ( rhs == 3 ) then
        v = varargin(3)
    end
    if ( rhs == 4 ) then
        m = varargin(3)
        n = varargin(4)
    end
    
    
    //
    // Check v,m,n
    //
    distfun_checkvmn ( "distfun_binornd" , 3 , varargin(3:$) )
    
    [N,Pr] = apifun_expandfromsize ( 2 , varargin(1:$) )
    if (N == []) then
        R = []
        return
    elseif(Pr == []) then
        R = []
        return
    end
    
    m = size(Pr,"r")
    n = size(Pr,"c")
    
    R = distfun_grandbin(m,n,N,Pr)
    
endfunction

distfun_binostat.sci

// Copyright (C) 2012 - Prateek Papriwal
//
// This file must be used under the terms of the CeCILL.
// This source file is licensed as described in the file COPYING, which
// you should have received as part of this distribution.  The terms
// are also available at
// http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt

function [M,V] = distfun_binostat(N,Pr)
    // Binomial  mean and variance
    //
    // Calling Sequence
    //   M = distfun_binostat(Pr)
    //   [M,V] = distfun_binostat(Pr)
    //
    // Parameters
    //   Pr : a 1x1 or nxm matrix of doubles, the probability of getting success in a Bernoulli trial 
    //   N : a 1x1 or nxm matrix of doubles , the total number of binomial trials . N belongs to the set {1,2,3,4,.......}
    //   M : a matrix of doubles, the mean
    //   V : a matrix of doubles, the variance
    //
    // Description
    //   Computes statistics from the Binomial distribution.
    //
    //   Any scalar input argument is expanded to a matrix of 
    //   doubles of the same size as the other input arguments.
    //
    // Examples
    //
    // // Check with expanded N
    // N = [10 100 1000 10000];
    // Pr = 0.1;
    // [m,v] = distfun_binostat(N,Pr);
    // me = [1 10 100 1000];
    // ve = [0.9 9 90 900];
    //
    // //Check with expanded Pr
    // N = 100;
    // Pr = [0.1 0.2 0.3 0.4];
    // [m,v] = distfun_binostat(N,Pr);
    // me = [10 20 30 40];
    // ve = [9 16 21 24];
    // 
    // //Check with both N and Pr expanded;
    // N = [10 100 1000 10000];
    // Pr = [0.1 0.2 0.3 0.4];
    // [m,v] = distfun_binostat(N,Pr);
    // me = [1 20 300 4000];
    // ve = [0.9 16 210 2400];
    //
    // Bibliography
    // http://en.wikipedia.org/wiki/Binomial_distribution
    //
    // Authors
    // Copyright (C) 2012 - Prateek Papriwal
    [lhs,rhs] = argn()
    apifun_checkrhs("distfun_binostat",rhs,2)
    apifun_checklhs("distfun_binostat",lhs,1:2)
    // Check type
    //
    apifun_checktype("distfun_binostat",N,"N",1,"constant")
    apifun_checktype("distfun_binostat",Pr,"Pr",2,"constant")
    
    //Check content
    // 
    apifun_checkgreq("distfun_binostat",N,"n",1,1)
    apifun_checkrange("distfun_binostat",Pr,"Pr",2,0,1)
    
    //
    [N,Pr] = apifun_expandvar(N,Pr)
    
    M = N .* Pr
    V = N .* Pr .* (1-Pr)
    
endfunction

2022-09-08 09:26