# 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))
}