# Extended binomial distribution

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

#########################################
#---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)
}

-
I have no knowledge of the language you are using but it seems to me that the first argument k of pebinom should be an integer and not a floating point number as in totalDensity = cumProb = pebinom(0.0,n,p,0.0). Also, it seems like you are including the point mass at $0$ twice in totalDensity and cumProb, once in the initialization and again when you call pebinom with k equal to $0$. –  Dilip Sarwate Nov 30 '11 at 15:37
Could you please supply a definition of the "extended binomial density"? Your function pebinom doesn't fit the bill because it does not define a valid probability distribution unless sum is identically 1. –  whuber Nov 30 '11 at 15:38
I have edited the code to comply with naming conventions. –  fbielejec Nov 30 '11 at 22:40

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.

-
We denote almost binomial distribution supported on $0,...,\left\lfloor N\right\rfloor +1$ by: $P\{X=k\}=\begin{cases} \left(\begin{array}{c} N\\ k \end{array}\right)p^{k}(1-p)^{N-k} & k=0,...,\left\lfloor N\right\rfloor \\ 1-\underset{i=1}{\overset{\left\lfloor N\right\rfloor }{\sum}}\left(\begin{array}{c} N\\ i \end{array}\right)p^{i}(1-p)^{N-i} & k=\left\lfloor N\right\rfloor +1 \end{cases}$ Now $EBin(N,p)$ has a density equal to $P\{X=k+1\}$ for $x\in[k,k+1)$ with point mass of $P(X=0)$ at $0$, thus is supported on $[0,...,\left\lfloor N\right\rfloor +1]$ –  fbielejec Nov 30 '11 at 17:11
@fbielejec Since $N$ is apparently not an integer, what is meant by $\binom{N}{k}$? Something to do with Gamma functions such as $\Gamma(N+1)/\Gamma(k+1)\Gamma(N-k+1)$? If so, how does your pebinom compute this since you seem to be using choose (n,k) which presumably has integer arguments and returns the binomial coefficient? –  Dilip Sarwate Nov 30 '11 at 17:22
choose(n,k) is defined for all real numbers n and integers k. For $k\geq1$ it is defined as $\frac{n\cdot(n-1)\cdot...\cdot(n-k+1)}{k!}$ –  fbielejec Nov 30 '11 at 17:46
@fbielejec: Ok, thanks for the definition. You random variable is definitely discrete, thus it cannot have "density equal to P{X=k+1} for x∈[k,k+1)" since your variable cannot have, for example, value of 1/2. The density is zero between integer points, all you have is just a bunch of atoms. Thus, the cdf is (as usual) just sum(pebinom(0:floor(x), N, p)) (note that you broken the naming conventions, your pebinom should really be debinom, it's a density wrt counting measure, and debinom should be pebinom). –  Anton Korobeynikov Nov 30 '11 at 20:29
@Anton Korobeynikov: I have added a new routines to my previous post, mostly following Your comments. If this works out as expected I will vectorize code for speed. –  fbielejec Dec 2 '11 at 9:01

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.

-