/*
**  PSUMBERN.SRC    Procedure that obtains the probability 
**                  distribution of a sum of independent
**                  but heterogeneous Bernoullis
**
**  by Victor Aguirregabiria
**  Last revision: February, 2003
**
**  FORMAT:
**      psum = psumbern(p0) ;
**
**  Inputs:
**
**      p0  - (n x k) matrix of probabilities 
**            Each column contains the probabilities
**            associated with n Bernouilli variables
**            Column j:
**              Pr(X1j=1) | Pr(X2j=1) | ... | Pr(Xnj=1)
**  Outputs:
**
**      psum - (n+1 x k) matrix with the probability
**             distributions of the sum of Bernoullis.
**             Collumn j contains the distribution of:
**                  Sj = X1j + X2j + ... + Xnj
**             psum[i,j] = Pr(Sj = i-1)
**
*/

proc (1) = psumbern(p0) ;
  local n, k, psum, nvar, xj, j, flag, conta ;
  n = rows(p0) ;
  k = cols(p0) ;
  if (n==1) ;
    psum = (1-p0) | p0 ;
  else ;
    psum = zeros(n+1,k) ;
    nvar = 2^n ;
    xj = zeros(n,1) ;
    j=2 ;
    do while (j<=nvar) ;
      psum[sumc(xj)+1,.] = psum[sumc(xj)+1,.]
         + prodc((p0.^xj).*((1-p0).^(1-xj)))' ; 
      if (xj[n]==0) ;
        xj = xj[1:n-1] | 1 ;
      else ;
        flag=0 ;
        conta=n-1 ;
        do while (flag==0)AND(conta>1) ;
          if (xj[conta]==0) ;
            xj = xj[1:1:conta-1] | 1 | zeros(n-conta,1) ;
            flag=1 ;
          endif ;
          conta=conta-1 ;
        endo ;
        if (conta==1)AND(flag==0) ;
          xj = 1 | zeros(n-conta,1) ;
        endif ;
      endif ;
      j=j+1 ;
    endo ;
    psum[sumc(xj)+1,.] = psum[sumc(xj)+1,.]
           + prodc((p0.^xj).*((1-p0).^(1-xj)))' ; 
  endif ;
  retp(psum) ;
endp ;


