# RProject5_HypothesisTesting.r

# 1.0   Two coins:  coin0 and coin1
#
#   P(Head | coin0)=0.5 and P(Head | coin1) =0.7
prob.coin0=0.5
prob.coin1=0.7

#help(sample)
# 2.0 Choose a coin at random, count number of heads in 10 tosses
prob.coin.random=sample(c(prob.coin0,prob.coin1), size=1)

x=rbinom(n=1, size=10, prob=prob.coin.random)

x
## [1] 6
# 3.  Likelihood Table

#     For each outcome of X (column 1) 
#       report Likelihood of coin0 (column2) and of coin1 (column 3)
# X= number of heads on 10 tosses of coin
list.x=0:10

pdf.coin0=dbinom(list.x, size=10,prob=prob.coin0)
pdf.coin1=dbinom(list.x, size=10,prob=prob.coin1)

likelihood.table=cbind(
  x=list.x, like.coin0=pdf.coin0, like.coin1=pdf.coin1,
  likeratio= pdf.coin0/pdf.coin1)

likelihood.table
##        x   like.coin0   like.coin1    likeratio
##  [1,]  0 0.0009765625 0.0000059049 165.38171688
##  [2,]  1 0.0097656250 0.0001377810  70.87787866
##  [3,]  2 0.0439453125 0.0014467005  30.37623371
##  [4,]  3 0.1171875000 0.0090016920  13.01838588
##  [5,]  4 0.2050781250 0.0367569090   5.57930823
##  [6,]  5 0.2460937500 0.1029193452   2.39113210
##  [7,]  6 0.2050781250 0.2001209490   1.02477090
##  [8,]  7 0.1171875000 0.2668279320   0.43918753
##  [9,]  8 0.0439453125 0.2334744405   0.18822323
## [10,]  9 0.0097656250 0.1210608210   0.08066710
## [11,] 10 0.0009765625 0.0282475249   0.03457161
# Plot pmf function of X under H0 and H1
par(mfcol=c(1,1))
plot(list.x, pdf.coin0, xlab="x", ylab="p(x | theta)", 
     ylim=c(0, max(c(pdf.coin0, pdf.coin1))))

points(list.x, pdf.coin1, col='red')
title(main="p(x | theta) for H0 (Black) and H1 (Red)")
title(sub=paste(
  "H0: theta = ", prob.coin0, ",  H1: theta = ", prob.coin1, collapse=""))

# Plot Likelihood Ratio H0:H1 as a function of X

par(mfcol=c(1,1))
plot(likelihood.table[,"x"], likelihood.table[,"likeratio"], xlab="x", ylab= "LikeRatio",
     main="Likelihood Ratio  H0:H1")
list.levels=c(1/10, 1,10)
abline(h=1, col='grey')
abline(h=list.levels, col=c(2,3,4))

plot(likelihood.table[,"x"], log(likelihood.table[,"likeratio"]), xlab="x", ylab= "Log LikeRatio",
     main="Likelihood Ratio  HO:H1")
abline(h=0,col='gray')
abline(h=log(list.levels), col=c(2,3,4))

paste(list.levels, collapse=", ")
## [1] "0.1, 1, 10"
title(sub=paste("LR Levels:  ", paste(list.levels, collapse=", "), sep=""))

# 3.0 Bayes Rule accepts H0 if Likelihood Ratio > c
#     where c=P(H1)/P(H0);   the prior odds of H1 to H0

# 4.0 Consider decision rule for each value of x*=0,1, ...,10
#
# d.x*=1 if x >=x*

#Create table of probability of errors
list.xstar=c(-1:10)

table.errorProbs<-matrix(NA, nrow=length(list.xstar), ncol=3)
for (j.xstar in c(1:length(list.xstar))){
  
  xstar=list.xstar[j.xstar]
  table.errorProbs[j.xstar,1]=xstar
  prob.rejectH0.given.H0=1-pbinom(xstar,size=10,prob=prob.coin0)
  table.errorProbs[j.xstar, 2]=prob.rejectH0.given.H0
  prob.acceptH0.givenH1 = pbinom(xstar,size=10,prob=prob.coin1)
  table.errorProbs[j.xstar, 3]=prob.acceptH0.givenH1
}
dimnames(table.errorProbs)[[2]]<-c("xstar", "P(Reject H0 | H0)", "P(Accept H0 | H1)")

print(table.errorProbs)
##       xstar P(Reject H0 | H0) P(Accept H0 | H1)
##  [1,]    -1      1.0000000000      0.0000000000
##  [2,]     0      0.9990234375      0.0000059049
##  [3,]     1      0.9892578125      0.0001436859
##  [4,]     2      0.9453125000      0.0015903864
##  [5,]     3      0.8281250000      0.0105920784
##  [6,]     4      0.6230468750      0.0473489874
##  [7,]     5      0.3769531250      0.1502683326
##  [8,]     6      0.1718750000      0.3503892816
##  [9,]     7      0.0546875000      0.6172172136
## [10,]     8      0.0107421875      0.8506916541
## [11,]     9      0.0009765625      0.9717524751
## [12,]    10      0.0000000000      1.0000000000
#


plot(table.errorProbs[,1], table.errorProbs[,2], 
     xlab="x* (Reject H0 if X > x*)",
    ylab="P(Reject H0 | H0)",
    main="P(Type I Error)")

par(mfcol=c(1,1))

plot(table.errorProbs[,1], table.errorProbs[,3], 
     xlab="x* (Reject H0 if X > x*)",
     ylab="P(Accept H0 |H1)",
     main="P(Type II Error)")

par(mfcol=c(1,1))

plot(table.errorProbs[,2], table.errorProbs[,3], 
     xlab="P(Reject H0 | H0)",
     ylab="P(Accept H0 | H1)",
     main="Risk Points of Decision Rules\nP(Type II Error) vs P(Type I Error)"
)

# Add Risk Points of decision rules d2.x* (opposite of rules d.x*)

# for constants c in list.xstar

# d.xstar=1 if x >=xstar
# d2.xstar=1 if x < xstar

table.errorProbs2<-matrix(NA, nrow=length(list.xstar), ncol=3)
for (j.xstar in c(1:length(list.xstar))){
  
  c=list.xstar[j.xstar]
  table.errorProbs2[,1]=c
  prob.rejectH0.given.H0=pbinom(c,size=10,prob=prob.coin0)
  table.errorProbs2[j.xstar, 2]=prob.rejectH0.given.H0
  prob.acceptH0.givenH1 = 1-pbinom(c,size=10,prob=prob.coin1)
  table.errorProbs2[j.xstar, 3]=prob.acceptH0.givenH1
}

points(table.errorProbs2[,2], table.errorProbs2[,3],pch="x",col="red")