Extended binomial distribution

Question

I'm trying to implement the extended binomial density function with support on c( 0 : (floor(N) + 1)), but I'm running into (I think) precision issues, as running:

########################
#---DENSITY FUNCTION---#
########################
debinom <- function(k, n, p, sum) {
    if (k <=  n) {
  return( choose(n, k) * p^k * (1-p)^(n-k) )
  } else {
    return (1.0 - sum)
    }
}#END: pebinom

##########################################
#---CUMULATIVE DISTRIBUTION FUNCTION 2---#
##########################################
pebinom <- function(x, n, p) {

  # point mass at 0
  totalDensity = cumProb = debinom(0.0, n, p, 0.0)

  k = 0
  while (k <= (x)) {
    density2 = debinom(k, n, p, totalDensity)
    totalDensity = totalDensity + density2
    cumProb = cumProb + density2
    k = k + 1
  }

  k = k + 1
  density = debinom(k, n, p, totalDensity)
  cumProb = cumProb + density * (x - k)

  return (cumProb) 
}#END: debinom

############
#---TEST---#
############
for (i in 0:10) {
x = i + runif(1)
cat(x, " ", pebinom(x, 100, 0.1), "\n")
}

gives a negative probabilities for tail values.

EDIT

I have changed, and mostly simplified the routines along the comments and answers I've received:

#########################################
#---PROBABILITY DISTRIBUTION FUNCTION---#
#########################################

debinom <- function(k, n, p) {

if (k <=floor(n)) {

  return( choose(n, k) * p^k * (1-p)^(n-k) )

  } else if(k == (floor(n)+1)) {

    cumProb = 0.0
    for(i in 0 : floor(n)) {
      cumProb = cumProb + debinom(i, n, p)  
    }

    return (1.0 - cumProb)

    } else {

  return(0.0)
    }

}#END: pebinom

########################################
#---CUMULATIVE DISTRIBUTION FUNCTION---#
########################################
pebinom <- function(x, N, P) {

cumProb = 0
for(i in 0 : (floor(x)) ) {
 cumProb = cumProb + debinom(i, N, P)
}

return(cumProb)
}

Answer 1

Well, it seems several problems here

First, what is "extended binomial distribution". Will you please give some reference to it? I'm aware about extended (aka truncated) negative binomial distribution, but here it seems it's not the same. Looking into pebinom() it looks like you're just looking for ordinary binomial distribution.

Second, the "else" case in pebinom looks completely bogus - if pebinom is a density function, when it should integrate (wrt counting measure) to one. Thus for k > n pebinom should return 0.

Next, couple of other problems: you're including the point mass at 0 twice, variables totalDensity and cumProb denotes the same thing.

Now, let's proceed to the problem of the negative values. It seems you're wanting to interpolate stuff if you x is between k and k + 1. In fact, you shouldn't do this, since the density wrt counting measure is zero between integer points. Anyway, your problem is that you're trying to interpolate at wrong points. Look, after your loop k is already > x, so, you have to interpolate between the points k-1 and k. Instead, you increase the k, calculate the density there and try to interpolate on a wrong interval, thus:

1. You should not increase the k after the loop (just remove that line)
2. You should multiply the density with (x - k + 1), not with x - k

PS: Get rid of loops, R perfectly allows you to write vectorized code.

Answer 2

There are several problems with your program. First, you set

totalDensity = cumProb = debinom(0.0, n, p, 0.0)

which sets totalDensity and cumProb to $(1-p)^n$. The first iteration of your loop (executed with k = 0) adds debinom(0, n, p, totalDensity) (which also has value $(1-p)^n$) to both totalDensity and cumProb and so now totalDensity and cumProb both contain $2(1-p)^n$. Thus, after the last iteration of the loop,

• totalDensity and cumProb have value $$(1-p)^n + \sum_{k=0}^{\lfloor x \rfloor} \binom{n}{k}p^k(1-p)^{n-k}$$ where the relationship, if any, between $x$ and $n$ is unspecified, but I am assuming that $x = n$. totalDensity and cumProb should have value equal to the sum but you have included an extra $(1-p)^n$. Thus, it is possible that totalDensity and cumProb might have value greater than $1$.

• k had value $\lfloor n \rfloor$ at the beginning of the last execution of the loop, and since k is incremented at the end of the loop, it has value $\lceil n \rceil$ upon exit from the loop.

You now increment k again to $\lceil n \rceil + 1$ and call debinom($\lceil$k$\rceil$, n, p, totalDensity) which returns 1 - totalDensity since the first argument of debinom is larger than the second, and this could well be a negative number. Then you multiply density by x-k which could well be a negative number if x is the same as n.

As a general comment, $$\binom{x}{k} = \frac{x(x-1)(x-2)\cdots (x-k+1)}{k!}$$ is positive for $0 \leq k \leq \lceil x \rceil$ and then alternates in sign as $k$ increases beyond $\lceil x \rceil$, and this might be a source of the negative probabilities being found at the tails.