#####################################################################################
#Brown Forsythe Test
#
#e is a vector of residuals, and x is some vector (usually yhat or one of the x_j's)

brown.forsythe=function(e,x){
  m=median(x)
  e1=e[x<=m]
  e2=e[x>m]
  e1med=median(e1)
  e2med=median(e2)
  d1=abs(e1-e1med)
  d2=abs(e2-e2med)
  d1bar=mean(d1)
  d2bar=mean(d2)
  n1=length(e1)
  n2=length(e2)
  sp=sqrt((sum((d1-d1bar)^2)+sum((d2-d2bar)^2))/(n1+n2-2))
  t=(d1bar-d2bar)/(sp*sqrt(1/n1+1/n2))
  return(2*pt(abs(t),n1+n2-2,lower.tail=FALSE))  
}


#######################################################################################
#Box Cox Transformation
#
#Y is the dependent variable from a regression equation, and X is the design matrix.
#lambdarange is a vector containing the possible values of lambda 

box.cox=function(Y,X,lambdarange){
  n=length(Y)
  K2=(prod(Y^(1/n)))
  L=length(lambdarange)
  SSE=1:L
  for(l in 1:L){
    lambda=lambdarange[l]
    K1=1/(lambda*K2^(lambda-1))
    if(lambda==0){W=K2*log(Y)}
    if(lambda!=0){W=K1*(Y^lambda-1)}
    tempmodel=lm(W~X)
    SSE[l]=deviance(tempmodel)
  }
  i=(sort(SSE,index.return=TRUE)$ix)[1]
  return(lambdarange[i])
}