{{ message }}

Instantly share code, notes, and snippets.

# coppeliaMLA/binDiff.R

Created Mar 21, 2014
A function that gives the probability mass function for the difference between to binomially distributed random variables
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 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 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 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 commented Nov 15, 2018

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