Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
A function that gives the probability mass function for the difference between to binomially distributed random variables
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))
@odelbor

This comment has been minimized.

Copy link

@odelbor odelbor commented Oct 27, 2016

Thanks for your helpful work. Please note a typo in the code: the two loops in diffBin should both start at i=0, rather than at i=1.

@markusloecher

This comment has been minimized.

Copy link

@markusloecher markusloecher commented Nov 14, 2018

Thanks for the code and also thanks @odelbor for pointing out the typo. I noticed that the distribution was not symmetric for the case p1=p2, n1=n2 and found the reason to be the wrong index i=1 start.

@markusloecher

This comment has been minimized.

Copy link

@markusloecher markusloecher commented Nov 15, 2018

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment