A function that gives the probability mass function for the difference between to binomially distributed random variables

binDiff.R

modBin<-function(k, n, p){
  
  if (k<=n) {
    return(dbinom(k, n, p))
  }
  else {
    return(0)
  }
}

diffBin<-function(z, n1, p1, n2, p2){

  prob<-0
  
  if (z>=0){  
    for (i in 1:n1){     
      prob<-prob+modBin(i+z, n1, p1)*modBin(i, n2, p2)
    }
    
  }
  else
  {
    for (i in 1:n2){     
      prob<-prob+modBin(i+z, n1, p1)*modBin(i, n2, p2)
    }
  }
  return(prob)
}

#Example
#n1=30, p1=0.5, n2=10, p2=0.7
n1<-20
p1<-0.1
n2<-10
p2<-0.7
s<-(-n2):n1

p<-sapply(s, function(x) diffBin(x, n1, p1, n2, p2))
plot(s,p)

#Montecarlo check

x<-rbinom(10000, n1, p1)
y<-rbinom(10000, n2, p2)
z<-x-y

t<-table(z)
lines(as.numeric(names(t)), t/sum(t))

Also, the dbinom function in R already (nicely) returns 0 outside its support.