#Discrete random variables and their distributions 1. Binomial dbinom(x=2,size=20,prob=0.5) pbinom(x=2,size=20,prob=0.5) qbinom(p=0.4,size=20,prob=0.5) rbinom(n=5,size=20,prob=0.5) plot(dbinom(0:20,size=20,prob=0.5),type="h") plot(dbinom(0:20,size=20,prob=0.8),type="h") 2. The Hypergeometric Distribution dhyper(x=2, m=10, n=30, k=6) phyper(q=2, m=10, n=30, k=6) qhyper(0.3, m=10, n=30, k=6) rhyper(nn=10, m=10, n=30, k=6) 3. The Geometric Distribution,let X count the number of failures before the first successes dgeom(4,prob=0.8) pgeom(4, prob = 0.8) qgeom(0.4,prob=0.8) rgeom(10,prob=0.8) plot(dgeom(0:20,prob=0.5),type="h") plot(dgeom(0:20,prob=0.8),type="h") 4.The Negative Binomial Distribution, let X count the number of failures before r successes dnbinom(x=5,size=3,prob=0.4) pnbinom(5,size=3,prob=0.4) qnbinom(0.5,size=3,prob=0.4) rnbinom(n=10,size=3,prob=0.4) plot(dnbinom(0:20,size=5,p=0.5),type="h") 5. Poisson distributino dpois(x=0,lambda=2.4) ppois(q=10,lambda=2.4) qpois(p=0.9,lambda=2.4) rpois(n=10,lambda=2.4) plot(dpois(0:20,lambda=1),type="h") x <- 0:20 plot(x, ppois(x, 1), type="s", lty=1,ylab="F(x)", main="Poisson approx of binomial") lines(x, pbinom(x, 100, 0.01),type="s",col=2,lty=2) legend("bottomright",legend=c("Poisson","Binomial"),lty=1:2,col=1:2) 6. Poisson and normal approximation of binomial probabilities, with estimated parameters #P(X<=k)=pbinom(k,n,p) #Poisson approximation: P(X<=k) app ppois(k,np) #Normal approximation: P(X<=k) app pnorm(k,np,npq) apprx <- function(n, p, R = 1000, k = 6) { trueval <- pbinom(k, n, p) # true binomial probability prob.zcc <- prob.zncc <- prob.pois <- NULL q<-1-p for (i in 1:R) { x <- rnorm(n, n * p, sqrt(n * p * q)) z.cc <- ((k + .5) - mean(x))/sd(x) # with cont. correction prob.zcc[i] <- pnorm(z.cc) z.ncc <- (k - mean(x))/sd(x) # no cont. correction prob.zncc[i] <- pnorm(z.ncc) y <- rpois(n, n * p) prob.pois[i] <- length(y[y <= k])/n } list(prob.zcc = prob.zcc, prob.zncc = prob.zncc, prob.pois = prob.pois, trueval = trueval) } R <- 1000 set.seed(10) out <- apprx(n = 200, p = .03, k = 6, R = 1000) # windows(6,5) plot(1:R, out$prob.pois, type = "l", col = "green", xlab = "Runs", main = expression(paste("Simulated Probabilities: ", n==200, ", ", p==0.03, sep="")), ylab = "Probability", ylim = c(.3, .7)) abline(h = out$trueval, col="red", lty=2) lines(1:R, out$prob.zcc, lty = 1, col = "purple") lines(1:R, out$prob.zncc, lty = 1, col = "orange") legend("bottomleft", c("Poisson", "Normal (with cc)", "Normal (w/o cc)"), lty = c(1), col = c("green", "purple", "orange")) set.seed(10) out <- apprx(n = 200, p = .03, k = 6, R = 1000) # windows(6,5) boxplot(out$prob.pois, boxwex = 0.25, at = 1:1 - .25, col = "green", main = expression(paste("Approximating Binomial Probability: ", n==200, ", ", p==0.03, sep="")), ylab = "Probablity", ylim = c(out$trueval - 0.2, out$trueval + 0.25)) boxplot(out$prob.zcc, boxwex = 0.25, at = 1:1 + 0, add = T, col = "purple") boxplot(out$prob.zncc, boxwex = 0.25, at = 1:1 + 0.25, add = T, col = "orange" ) abline(h = out$trueval, col = "red", lty=2) legend("topleft", c("Poisson", "Normal (with cc)", "Normal (w/o cc)"), fill = c("green", "purple", "orange"))