linear_regression_indicator_matrix = function(XTrain,yTrain){
  #
  # R code to compute classification based on linear regression.
  #
  # See the section with the same name as this function in Chapter 4 from the book ESLII
  # Note one could prob use the R command "lm" to do this routine 
  #
  # Inputs:
  #   XTrain = training matrix (without the common column of ones needed to represent the constant offset)
  #   yTrain = training labels matrix of true classifications with indices 1 - K (where K is the number of classes)
  #
  # Written by:
  # -- 
  # John L. Weatherwax                2009-04-21
  # 
  # email: wax@alum.mit.edu
  # 
  # Please send comments and especially bug reports to the
  # above email address.
  # 
  #-----

  K = max(yTrain) # the number of classes
  N = dim( XTrain )[1] # the number of samples 

  # form the indicator responce matrix Y
  # 
  Y = mat.or.vec( N, K ) 
  for( ii in 1:N ){
    jj       = yTrain[ii] 
    Y[ii,jj] = 1.
  }
  Y = as.matrix(Y)

  # for the feature vector matrix XTrain ( append a leading column of ones ) and compute Yhat:
  #
  ones = as.matrix( mat.or.vec( N, 1 ) + 1.0 )
  Xm   = as.matrix( cbind( ones, XTrain ) )
  Bhat = solve( t(Xm) %*% Xm, t(Xm) %*% Y ) # this is used for predictions on out of sample data  
  Yhat = Xm %*% Bhat # the discriminant predictions on the training data  
  gHat = apply( Yhat, 1, 'which.max' ) # classify this data 
  
  return( list(Bhat,Yhat,gHat) )

}