#### verify Ex7.7.7(b)
## verify the density function 
n <- 10                           # sample size
xbar <- 1                         # sample mean
S <- 2                            # sample standard deviation
x1 <- xbar - (n-1)*S/sqrt(n)      # left end of x1 given xbar and S
x2 <- xbar + (n-1)*S/sqrt(n)      # right end of x1 given xbar and S

xf <- function(x) {               # density function of x1 given xbar and S
  gamma((n-1)/2)/(gamma(1/2)*gamma(n/2-1))*sqrt(n)/((n-1)*S)*(1-((x-xbar)/((n-1)*S/sqrt(n)))^2)^(n/2-2);
}

integrate(xf, lower=x1, upper=x2) # it should give us 1 as the total probability

## verify the MVUE
n <- 10                           # sample size
theta1 <- 1                       # mean of normal population
theta2 <- 4                       # variance of normal population
X <- rnorm(n, mean=theta1, sd=sqrt(theta2))  # simulate a random sample form N(theta1, theta2)
xbar <- mean(X)                   # sample mean
S <- sqrt(var(X))                 # sample standard deviation

fmle <- function(c) {             # mle for P(X<c)
  pnorm((c-xbar)/S)
}

fmvue <- function(c) {            # MVUE for P(X<c)
  if(c <= (xbar-(n-1)*S/sqrt(n))) ans <- 0 else {
    if(c <= xbar) ans <- 0.5 - 0.5*pbeta(((c-xbar)/((n-1)*S/sqrt(n)))^2, shape1=0.5, shape2=n/2-1) else {
      if(c <= (xbar+(n-1)*S/sqrt(n))) ans <- 0.5 + 0.5*pbeta(((c-xbar)/((n-1)*S/sqrt(n)))^2, shape1=0.5, shape2=n/2-1) else {
        ans <- 1
      }
    }
  }
  ans
}

tvec <- seq(-10,10,by=0.01)       # c values from -10 to 10
m <- length(tvec)                 # number of different c values
mlevec <- mvuevec <- rep(0, m)    
for(i in 1:m) {
  mlevec[i] <- fmle(tvec[i]);     # mle for P(X < tvec[i])
  mvuevec[i] <- fmvue(tvec[i]);   # MVUE for P(X < tvec[i])
}

par(mfrow=c(1,1))                 # plot mle and MVUE for different c values
plot(tvec,mlevec,type="l",col="blue",lty=1,lwd=2,xlab="c",ylab="mle or MVUE", main="Estimates of P(X<c)")
points(tvec,mvuevec,type="l",col="red",lty=1,lwd=2)
legend(tvec[1],max(mlevec),lty=c(1,1),lwd=c(2,2),col=c("blue","red"),legend=c("mle","MVUE"))