# Binomial -------------------------------------------------------------------------------

# What is the chance of getting exactly 2 sixes on 10 rolls of a die?

my.binomial <- function(k,n,p){
(choose(n,k)* p^k * (1-p)^(n-k))
}

my.binomial(2,10,1/6)

# or R already has a function for this
dbinom(2,10,1/6)

# Exploring Binomial Distribution---------------------------------------------------------

dice <- function(nroll=60) {
sum(runif(nroll,0,1) < 1/6)
}

diceout <- function(nsim=1000){
dice1 <<- rep(NA,nsim)
for (i in 1:nsim)
{dice1[i] <<- dice()}
hist(dice1,breaks=0:60)
}

# Poisson --------------------------------------------------------------------------------

# What is the chance of seeing 18 cases of gonorrhea in a county when we expect to see 7?

my.poisson  <- function(l,k){
l^k*exp(-l)/factorial(k)
}

my.poisson(7,18)

# or R already has a function for this too
dpois(18,7)

# density for each value greater than 17....
dpois(18:25,7)
# can calculate tail probabilities based on this
sum(dpois(18:50,7))

# or R has a function for this too
ppois(17,7,lower.tail=FALSE)

# Comparing the two ----------------------------------------------------------------------

n <- 10000
p <- .0005

# fine when n is large and p is small

dbinom(2,n,p)
dpois(2,n*p)
dpois(2,5)

# or tail probabilities

pbinom(2,n,p)
ppois(2,n*p)

# but not probability of getting exactly no heads on 3 flips of a coin

n <- 3
p <- 1/2

dbinom(0,n,p)
dpois(0,n*p)

# Use Poisson Regression to calculate rate ratios ------------------------------------------------

pop          <- c(2397329,16095710)
cases        <- c(6643,4877)
race         <- c(1,0)
names(race)  <- c("Black","White")
dat          <- data.frame(cases,pop,race)

rate <- cases/pop
rr   <- rate[1]/rate[2]

p.out <- glm(cases~race,family=poisson(link="log"),offset=log(pop))
rr.p  <- exp(p.out$coef[2])

rr
rr.p