# Section 4; Question 1:
#
n = 86
mu0 = 494
s = 124
ybar = 502
mu = 500
alpha = 0.05
z_alpha_over_two = qnorm( 1-alpha/2 )
beta = pnorm( (mu0 - mu)/(s/sqrt(n)) - z_alpha_over_two ) + (1-pnorm( (mu0 - mu)/(s/sqrt(n)) + z_alpha_over_two ))
print(sprintf('beta= %f', beta))


# Section 4; Question 2:
#
n = 30
mu0 = 25
s = 2.4
alpha = 0.1
z_alpha = qnorm( 1 - alpha )
print( mu0 + (s/sqrt(n)) * z_alpha )


# Section 4; Question 3:
#
z_crit = qnorm(1-0.06/2) # taken from Section 2; Question 2:
y_limits = 95 + (15/sqrt(22)) * z_crit * c(-1, +1)
print(y_limits)

mu = 90
beta = pnorm( (y_limits[2] - mu)/(15/sqrt(22)) ) - pnorm( (y_limits[1] - mu)/(15/sqrt(22)) )
print(sprintf('1-beta= %f', 1-beta))


# Section 4; Question 4:
#
alpha = 0.05
mu0 = 60
n = 16
s = 4
z_alpha_over_two = qnorm( 1-alpha/2 )
mus = seq(50, 70, length.out=100 )
beta = pnorm( (mu0 - mus)/(s/sqrt(n)) + z_alpha_over_two ) - pnorm( (mu0 - mus)/(s/sqrt(n)) - z_alpha_over_two )
#postscript("../../WriteUp/Graphics/Chapter6/chap_6_4_4_power_plot.eps", onefile=FALSE, horizontal=FALSE)
plot( mus, 1-beta, type='l', xlab='mu', ylab='power=1-beta' )
grid()
#dev.off()


# Section 4; Question 5:
#
mu0 = 240
alpha = 0.01
n = 25
s = 50
mu = 220
z_alpha = qnorm( 1-alpha )
beta = 1 - pnorm( (mu0-mu)/(s/sqrt(n)) - z_alpha )
print(sprintf('beta= %f', beta))


# Section 4; Question 6:
#
n = 36
s = 8.0
mu0 = 60
alpha = 0.07
y_bar_star = (s/sqrt(n)) * qnorm((alpha+1)/2)
print(sprintf('y_bar_star= %f', y_bar_star))

mu = 62
power = pnorm( (60+y_bar_star-mu)/(s/sqrt(n)) ) - pnorm( (60-y_bar_star-mu)/(s/sqrt(n)) )
print(sprintf('power (initial test): %f', power))

z_crit = qnorm(1-alpha/2)
y_limits = mu0 + (s/sqrt(n)) * z_crit * c(-1, +1)
print(y_limits)

power = pnorm( (y_limits[2]-mu)/(s/sqrt(n)) ) - pnorm( (y_limits[1]-mu)/(s/sqrt(n)) )
print(sprintf('power (correct test): %f', power))


# Section 4; Question 7:
#
mu0 = 200
alpha = 0.1
z_alpha = qnorm( 1-alpha )
s = 15.0
power = 0.75
z_power = qnorm( power )
mu = 197
sn = ( (z_power + z_alpha) * s )/ ( mu0 - mu )
print(sprintf('z_power= %f; n= %f', z_power, sn^2))


# Section 4; Question 8:
#
n = 45
mu0 = 10
alpha = 0.05
mu = 12
z_alpha_over_two = qnorm( 1 - alpha/2 )
s = 4
beta = pnorm( (mu0 - mu)/(s/sqrt(n)) + z_alpha_over_two ) - pnorm( (mu0 - mu)/(s/sqrt(n)) - z_alpha_over_two )
print(sprintf('n= %d; beta= %f', n, beta))


# Section 4; Question 9:
#
mu0 = 30
n = 16
power = 0.85
mu = 34
s = 9
z_beta = qnorm( 1-power )
z_alpha = z_beta - (mu0 - mu)/(s/sqrt(n))
alpha = qnorm(z_alpha)
print(sprintf('z_beta= %f; z_alpha= %f; alpha= %f', z_beta, z_alpha, alpha))


# Section 4; Question 10:
#
print(exp(-3.2))

print(1-exp(-0.75*3.2))

#postscript("../../WriteUp/Graphics/Chapter6/chap_6_4_10_pdf_plots.eps", onefile=FALSE, horizontal=FALSE)
ys = seq(0.01, 5, length.out=500)
pdf_1 = exp(-ys)
lam = 4/3
pdf_2 = (1/lam) * exp(-ys/lam)
plot( ys, pdf_1, type='l', xlab='y', ylab='f_Y(y)', col='black', xlim=c(0, 4), lwd=2 )
lines( ys, pdf_2, type='l', col='darkblue', lwd=2 )
polygon( c( ys[ys>3.2], 3.2 ), c( pdf_1[ys>3.2], 0 ), col=rgb(0, 0, 0, 0.5) )
polygon( c( 0, ys[ys<3.2], 3.2 ), c( 0, pdf_2[ys<3.2], 0 ), col=rgb(0, 0, 1, 0.25) )
abline(v=3.2, col='red')
legend('top', c('H_0', 'H_1'), col=c('black', 'blue'), lty=c(1, 1))
grid()
#dev.off()


# Section 4; Question 11:
#
alpha = 9/(131+9)
beta = 15/(125+15)
print(sprintf('alpha= %f; beta= %f', alpha, beta))


# Section 4; Question 12:
#
ks = 2:3
print(sprintf('alpha = %f', sum(dhyper(ks, 5, 5, 3))))
print(sprintf('power(w=6,r=4)= %f', sum(dhyper(ks, 6, 4, 3)) ))
print(sprintf('power(w=7,r=3)= %f', sum(dhyper(ks, 7, 3, 3)) ))


# Section 4; Question 13:
#
k = ( 2^5 * ( 1 - 0.05 ) )^(1/5)
print(sprintf('k= %f', k))


# Section 4; Question 16:
#
p_5 = dbinom( 1, 2, 1/2 ) * dbinom( 4, 4, 1/2 ) + dbinom( 2, 2, 1/2 ) * dbinom( 3, 4, 1/2 )
p_6 = dbinom( 2, 2, 1/2 ) * dbinom( 4, 4, 1/2 )
print(p_5)
print(p_6)
print(p_5 + p_6)


# Section 4; Question 18:
#
alpha = sum(dpois( 0:2, 6 ))
print(alpha)
beta = 1- sum(dpois( 0:2, 4 ))
print(beta)


# Section 4; Question 19:
#
ks = 1:3
p = 1/2
beta = sum( (1-p)^(ks-1) * p )
print(beta)


# Section 4; Question 22:
#
lhs_fn = function(x){
    theta = 2
    alpha = 0.05
    (1/theta^2) * ( 1 - log(x/theta^2) ) * x - alpha
}

# Plot the function we want to find the zero of:
#
ks = seq( 0.001, 0.15, length.out=100 )
lhs = lhs_fn(ks)
plot( ks, lhs, type='l', xlab='k star', ylab='left-hand-side' ); grid()
abline(h=0, col='green')

library(stats)
res = uniroot( lhs_fn, c(0.001, 0.05) )
print(sprintf('k_star= %f', res$root))