bstrap.nonparam = function(x, B, s, ...)
# T = bstrap.nonparam(x, B, s, ...)
#
# Nonparametric bootstrap: generate B bootstrap samples and for
# each compute the bootstrap replicate of the parameter of
# interest s(x, ...).
# Return the estimates in a vector T with B elements.
{
  n = length(x)
  bsamples = matrix(sample(x, size=n*B, replace=TRUE), nrow=B)
  theta.star = apply(bsamples, 1, s, ...)

  return (theta.star)
}

bstrap.theta0 = function(T)
# theta0 = bstrap.theta0(T)
#
# Estimate the mean value of the boostrap replicates of a
# parameter.
{
  return (mean(T))
}

bstrap.se = function(T)
# se = bstrap.se(T)
# 
# Standard error of a parameter whose bootstrap replicates are
# given in T.
{
  return (sd(T))
}

bstrap.bca = function(x, B, s, ..., alpha=c(0.025, 0.05)) 
# ci.bounds = bstrap.bca(x, B, s, ..., alpha)
#
# Bias-corrected and accelerated bootstrap (2*alpha)-level 
# confidence intervals for a parameter s(x,...), obtained 
# from B bootstrap replicates.
{
  n = length(x)
  
  theta.star = bstrap.nonparam(x, B, s, ...)
  
  theta.hat = s(x, ...)

  z0 = qnorm( sum(theta.star < theta.hat) / B )

  theta.hat.i = rep(0, n)
  for (i in 1:n) theta.hat.i[i] = s(x[-i], ...)

  d = mean(theta.hat.i) - theta.hat.i
  a.hat = sum(d^3) / (6 * (sum(d^2))^1.5)

  z.alpha.lo = qnorm(alpha)
  z.alpha.hi = qnorm(1-alpha)
  
  alpha.1 = pnorm(z0 + (z0 + z.alpha.lo)/(1 - a.hat * (z0 + z.alpha.lo)))
  alpha.2 = pnorm(z0 + (z0 + z.alpha.hi)/(1 - a.hat * (z0 + z.alpha.hi)))
  
  ci.lo = quantile(theta.star, probs=alpha.1)
  ci.hi = quantile(theta.star, probs=alpha.2)
  
  # alternative:
  # tt = sort(theta.star)
  # ci.lo = tt[c(trunc(B*alpha.1)]
  # ci.hi = tt[c(trunc(B*alpha.2)]

  return (list(ci.lo=ci.lo, ci.hi=ci.hi))
}

bstrap.tst = function(z, y, B, tst)
# asl.boot = bstrap.tst(z, y, B, tst)
#
# Bootstrap hypothesis testing procedure: use the test statistic tst()
# to compare the distributions of z and y. Returns the bootstrap estimate
# of ASL (achieved significance level) after B replicates.
{
  n = length(z)
  m = length(y)
  
  tst.0 = tst(z, y)
  
  x = c(z,y)
  bsamples = matrix(sample(x, size=(n+m)*B, replace=TRUE), nrow=B)
  tst.star = apply(bsamples, 1, function(x_) (tst(x_[1:n],x_[(n+1):(n+m)])) )
  
  return (sum(tst.star >= tst.0)/B)
}

tst.student = function(z, y)
# Implemented to be similar with the slides
{
  n = length(z)
  m = length(y)

  mu.z = mean(z)
  mu.y = mean(y)

  sigma.bar = sqrt(((n-1)*var(z) + (m-1)*var(y))/(m+n-2))
  
  return ( (mu.z - mu.y) / (sigma.bar*sqrt(1/n+1/m)) )
}

tst.meandiff = function(z, y)
# Implemented to be similar with the slides
{
  mu.z = mean(z)
  mu.y = mean(y)
  
  return ( mu.z - mu.y )
}