# BOUNDED LOGISTIC REGRESSION ESTIMATION USING A QUADRATIC PENALTY.  Returns
# the maximum penalized likelihood estimate for the parameters of a logistic
# regression model with bounded probabilities.  The inputs in the training 
# cases are in X, which should be a matrix (or convertable to a matrix) with 
# one row per training case, and one column per input variable.  The training 
# targets are in y, which should be a vector of 0/1 values.  The coefficients 
# (but not the intercept) are penalized by lambda (default 0) times the sum of 
# their squares.  The usual logistic probabilities are shifted and scaled so
# that there is a minimum and maximum probability for class 1, whose values are
# estimated as additional model parameters (in log form). 
#
# The first two elements of the parameter vector returned are the logs of the
# minimum probabilities for class 0 and for class 1; the third element is the 
# intercept; the remaining elements are the coefficients of the inputs.  Initial
# values for the parameters may be specified in the same order (default is -1
# for the log minimum probabilities and 0 for other parameters).  The initial
# value is extended with zeros if shorter than the total number of parameters,
# and truncated to the needed length if it is longer than needed.

lrb.est = function (y, X, lambda, init=c(-1,-1))
{
  X = as.matrix(X)
  init = c(init, rep(0, ncol(X)+3)) [1:(ncol(X)+3)]

  penalized.minus.log.like = 
    function (params) - lrb.log.like(y,X,params) + lambda*sum(params[-(1:3)]^2)

  result = nlm (penalized.minus.log.like, init, iterlim=200)

  params = result$estimate
  attr(params,"max.log.lik") = -result$minimum  # tack on max value as attr

  params
}

# BOUNDED LOGISTIC REGRESSION LOG LIKELIHOOD.  Computes the log probability of 
# the training data with inputs X and targets y, using the parameters in params.
# See the documentation on lrb.est for the format of the parameter vector.

lrb.log.like = function (y, X, params)
{
  X = as.matrix(X)
  ys = 2*y-1

  lin = as.vector(params[3] + X %*% params[-(1:3)])

  m = 0.5/(1+exp(-params[1:2]))

  sum (log (m[1+y] + (1-sum(m))/(1+exp(-lin*ys))))
}


# MAKE PREDICTIONS USING BOUNDED LOGISTIC REGRESSION.  Returns the probability 
# that y=1 for each of the test cases whose inputs are in X, using the 
# parameters in params.  See the documentation on lrb.est for the format of 
# the parameter vector.

lrb.pred = function (X, params)
{
  X = as.matrix(X)

  lin = as.vector(params[3] + X %*% params[-(1:3)])

  m = 0.5/(1+exp(-params[1:2]))

  m[2] + (1-sum(m))/(1+exp(-lin))
}