# pset 7

#problem 1 ----
#1a ----
n = 250
x = 140
theta.0 = .5
p1 = 1-pbinom(x-1,n,theta.0)
s = sprintf("one-sided p = %.5f", p1)
print(s)
s = sprintf("one-sided p = %.5f", 2*p1)
print(s)

#1b ----
n = 250
alpha1 = .05
theta.0 = .5
#cp stands for critical point
cp1.left = qbinom(alpha1/2,n,theta.0) -1
cp1.right = qbinom(1- alpha1/2,n,theta.0) + 1
print(cp1.left)
print(cp1.right)
#Check
sum(dbinom(0:cp1.left,n,theta.0))
sum(dbinom(cp1.right:n,n,theta.0))
# widen by 1 and see the probability > .005 in each side
sum(dbinom(0:(cp1.left+1),n,theta.0))
sum(dbinom((cp1.right-1):n,n,theta.0))

alpha2 = .1
#cp stands for critical point
cp2.left = qbinom(alpha2/2,n,theta.0) -1
cp2.right = qbinom(1- alpha2/2,n,theta.0) + 1
print(cp2.left)
print(cp2.right)
#Check
sum(dbinom(0:cp2.left,n,theta.0))
sum(dbinom(cp2.right:n,n,theta.0))
# widen by 1 and see the probability > .005 in each side
sum(dbinom(0:(cp2.left+1),n,theta.0))
sum(dbinom((cp2.right-1):n,n,theta.0))

# make a plot # NEEDS n and theta.0 already defined
x = 80:180
prob.0 = dbinom(x,n,theta.0)
par('mar'=c(5.1,4.1,1.1,1.1))
plot(x, prob.0, type='p', pch=19, col='blue',ylab="prob",ylim=c(-.01,.06), cex=.8, cex.lab=1.2)
col1 = rgb(1,.65,.30)
col2 = rgb(0,.64,.84)
lines(c(80,cp1.left),c(-.005,-.005), col=col1, lwd=8)
lines(c(cp1.right,180),c(-.005,-.005), col=col1, lwd=8)
lines(c(80,cp2.left),c(-.008,-.008), col=col2, lwd=8)
lines(c(cp2.right,180),c(-.008,-.008), col=col2, lwd=8)
text(150,.05,  bquote(theta == .(theta.0)), cex=1.6)

#1c ----
n = 250
alpha = .01
theta.0 = .5
criticalPoint.left = qbinom(alpha/2,n,theta.0) -1
criticalPoint.right = qbinom(1- alpha/2,n,theta.0) + 1
print(criticalPoint.left)
print(criticalPoint.right)
#Check
sum(dbinom(0:criticalPoint.left,n,theta.0))
sum(dbinom(criticalPoint.right:n,n,theta.0))
# widen by 1 and see the probability > .005 in each side
sum(dbinom(0:(criticalPoint.left+1),n,theta.0))
sum(dbinom((criticalPoint.right-1):n,n,theta.0))

#1d ----
n = 250
alpha = .05
theta.0 = .5
criticalPoint.left = qbinom(alpha/2,n,theta.0) -1
criticalPoint.right = qbinom(1- alpha/2,n,theta.0) + 1
print(criticalPoint.left)
print(criticalPoint.right)
#Check
sum(dbinom(0:criticalPoint.left,n,theta.0))
sum(dbinom(criticalPoint.right:n,n,theta.0))
# widen by 1 and see the probability > .005 in each side
sum(dbinom(0:(criticalPoint.left+1),n,theta.0))
sum(dbinom((criticalPoint.right-1):n,n,theta.0))
# Compute power
rejectRegion = c(0:criticalPoint.left, criticalPoint.right:n)
theta.55 = .55
power.55 = sum(dbinom(rejectRegion,n,theta.55))
print(power.55)
theta.6 = .6
power.6 = sum(dbinom(rejectRegion,n,theta.6))
print(power.6)

#d.ii
x = 80:180
n = 250
theta.0 = .5
prob.0 = dbinom(x,n,theta.0)

theta = .55
prob.A = dbinom(x,n,theta)
par('mar'=c(5.1,4.1,1.1,1.1))
plot(x, prob.0, type='p', pch=19, col='blue',ylab="prob",ylim=c(-.01,.06),
     cex=.8, cex.lab=1.2)
points(x, prob.A, type='p', pch=19, col='red',ylab="prob", cex=.8)
lines(c(80,109),c(-.005,-.005), col='green', lwd=8)
lines(c(141,180),c(-.005,-.005), col='green', lwd=8)
text(170,.05,  bquote(theta == .(theta)), cex=1.6)

theta = .6
prob.A = dbinom(x,n,theta)
par('mar'=c(5.1,4.1,1.1,1.1))
plot(x, prob.0, type='p', pch=19, col='blue',ylab="prob",ylim=c(-.01,.06),
     cex=.8, cex.lab=2)
points(x, prob.A, type='p', pch=19, col='red',ylab="prob", cex=.8)
lines(c(80,109),c(-.005,-.005), col='green', lwd=8)
lines(c(141,180),c(-.005,-.005), col='green', lwd=8)
text(170,.05,  bquote(theta == .(theta)), cex=1.6)

#1e ----
theta = .55
# Here's code to do this for one value of n --see below for loop
n = 1055
criticalPoint.left = qbinom(.025,n,.5) - 1;
criticalPoint.right = qbinom(.975,n,.5) + 1;
cat(criticalPoint.left, '  ', criticalPoint.right)
rejectionRegion = c(0:criticalPoint.left, criticalPoint.right:n)
power = sum(dbinom(rejectionRegion, n, theta))
print(power)

#Here we search through values of n until we get power of .9
#We do a naive search because the power is not strictly
#monotonic in n. This is because the exact significance
#is always a small amount less than the nominal
#significance of .05. The amount less can fluctuate
#with n. Likewise the power can fluctuate a little.
theta = .55
n = 250
done = FALSE;
while (!done)
{
  #find critical points for rejection region (based on theta=.5)
  criticalPoint.left = qbinom(.025,n,.5) - 1;
  criticalPoint.right = qbinom(.975,n,.5) + 1;
  rejectionRegion = c(0:criticalPoint.left, criticalPoint.right:n)
  power = sum(dbinom(rejectionRegion, n, theta))
  cat(n,'  ',power,'\n')

  #Check for end condition
  if (power >= .9)
    done = TRUE
  else
    n = n+1
}

# Plot for 1e
x = 470:630
n = 1055
theta.0 = .5
alpha = .05
prob.0 = dbinom(x,n,theta.0)

theta = .55
prob.A = dbinom(x,n,theta)
par('mar'=c(5.1,4.1,1.1,1.1))
plot(x, prob.0, type='p', pch=19, col='blue',ylab="prob",ylim=c(-.005,.03),
     cex=.6, cex.lab=1.2)
points(x, prob.A, type='p', pch=19, col='red',ylab="prob", cex=.6)
lines(c(470,criticalPoint.left),c(-.005,-.005), col='green', lwd=8)
lines(c(criticalPoint.right,630),c(-.005,-.005), col='green', lwd=8)
text(500,.028,  bquote(theta == .(theta)), cex=1.6)
text(560,.028,  bquote(n == .(n)), cex=1.6)

#1f ----
hyp = c(.5,.55)
prior = c(.5,.5)
n = 250
data = 140
lik = dbinom(data,n,hyp)
print(lik)
posterior.unnorm = prior*lik
posterior = posterior.unnorm/sum(posterior.unnorm)
print(posterior)

#2a ----
p.typeI = 9/140
p.typeII = 15/140
print(p.typeI)
print(p.typeII)

# 3a ---
m = 11
s2 = 4
m0 = 10
n = 16
t = (m-m0)/sqrt(s2/n)
print(t)
p = 2*pt(-abs(t), n-1)
print(p)

#3b
print(p/2)

#4b----
thresh = qnorm(.96,40,5/sqrt(3))
print(thresh)

#4c ----
n = 3
mu = 40
sigma = 5/sqrt(n)
alpha = .04
xcrit = qnorm(1-alpha, mu, sigma)
power = 1-pnorm(xcrit,45,sigma)
print(power)

# Compute power for various n
for (n in 4:16)
{
  sigma = 5/sqrt(n)
  alpha = .04
  xcrit = qnorm(1-alpha, mu, sigma)
  power = 1-pnorm(xcrit,45,sigma)
  s = sprintf("%d, %f, %f", n, power, xcrit)
  print(s)
}