# SET THE WORKING DIRECTORY
setwd("C:/Courses/STA5505/Summer2008/")



# PLOT OF POISSON PROBABILITIES
mu <- c(1, 5, 10, 100)
par(mfrow=c(2,2), mar=c(4,4,4,4))
for (i in 1:4){
    x <- (max(0, mu[i]-50) : (mu[i]+50))
    y <- dpois(x, lambda=mu[i]) 
    plot(x, y, ylab="probablity", xlab="x", 
        main=paste("mu=", mu[i], sep=""), type="p")
}



library(combinat)
 n <- c(100,20,10)
 p <- matrix(c(.3,.1,.5,.1,.1,.2,.6,.8,.3),3)
 my.rmultinomial(n,p)
 
my.rmultinomial<-function (n, p, rows = max(c(length(n), nrow(p)))) 
{
    rmultinomial.1 <- function(n, p) {
        k <- length(p)
        tabulate(sample(k, n, replace = TRUE, prob = p), nbins = k)
    }
    n <- rep(n, length = rows)
    p <- p[rep(1:nrow(p), length = rows), , drop = FALSE]
    sapply(1:rows, function(i) rmultinomial.1(n[i],  p[i, ]))
}


##### Vegetarian Example ######

# R Code for Wilson's confidence interval for a proportion (Journal of the American Statistical 
# Association 1927). This is the score CI, based on inverting the asymptotic normal test 
# using the null standard error. 

scoreci <- function(x,n,conflev)
{
  zalpha <- abs(qnorm((1-conflev)/2))
  phat <- x/n
  bound <- (zalpha*((phat*(1-phat)+(zalpha**2)/(4*n))/n)**(1/2))/
           (1+(zalpha**2)/n)
  midpnt <- (phat+(zalpha**2)/(2*n))/(1+(zalpha**2)/n)

  uplim <- round(midpnt + bound,digits=4)
  lowlim <- round(midpnt - bound,digits=4)

  results <- data.frame(lowlim,uplim)

  cat("\n")
  cat("With confidence level",conflev," and sample proportion",
     round(phat,digits=4),
     " \nthe lower and upper limits for the score confidence interval are: \n")
  cat("\n")
  print(results)
  cat("\n")

  # This function computes a confidence interval for a proportion.  
  # It is based on inverting the large-sample normal score test for the
  # proportion.  The input parameters are the number of successes, x, the
  # total sample size, n, and the confidence level, conflev (e.g. 0.90).  
  # Returned are the endpoints for the 100*conflev % confidence
  # interval for the proportion.

  # binconf(x,n,conflev)
}
scoreci(x=0,n=25,conflev=0.95)

# prop.test
res<-prop.test(x=0,n=25,correct=F)
res$conf.int

## exact intervals
binom.test(x=0, n=25)



######################################################
### Pearson's Chi-square Statistic (Mendel's theories)
chisq.test(x=c(6022,2001),p=c(.75,.25))

## LR chi-squared test

# mendel's theories
obs<-c(6022,2001)
expected<-8023*c(.75,.25)
1-pchisq(2 * sum(obs * log(obs/expected)),df=1)