\HeaderA{pcls}{Penalized Constrained Least Squares Fitting}{pcls}
\keyword{models}{pcls}
\keyword{smooth}{pcls}
\keyword{regression}{pcls}
\begin{Description}\relax
Solves least squares problems with quadratic penalties subject to linear
equality and inequality constraints using quadratic programming.
\end{Description}
\begin{Usage}
\begin{verbatim}
pcls(M)
\end{verbatim}
\end{Usage}
\begin{Arguments}
\begin{ldescription}
\item[\code{M}] is the single list argument to \code{pcls}. It should have  the 
following elements:
\Itemize{
\item[y] The response data vector.
\item[w] A vector of weights for the data (often proportional to the 
reciprocal of the variance). 
\item[X] The design matrix for the problem, note that \code{ncol(M\$X)}
must give the number of model parameters, while \code{nrow(M\$X)} 
should give the number of data.
\item[C] Matrix containing any linear equality constraints 
on the problem (e.g. \eqn{ \bf C}{C} in \eqn{ {\bf Cp}={\bf
c} }{Cp=c}). If you have no equality constraints
initialize this to a zero by zero matrix. Note that there is no need 
to supply the vector \eqn{ \bf c}{c}, it is defined implicitly by the 
initial parameter estimates \eqn{ \bf p}{p}.
\item[S] A list of penalty matrices. \code{S[[i]]} is the smallest contiguous matrix including 
all the non-zero elements of the ith penalty matrix. The first parameter it
penalizes is given by \code{off[i]+1} (starting counting at 1). 
\item[off] Offset values locating the elements of \code{M\$S} in
the correct location within each penalty coefficient matrix. (Zero
offset implies starting in first location)
\item[sp] An array of  smoothing parameter estimates.
\item[p] An array of feasible initial parameter estimates - these must
satisfy the constraints, but should avoid satisfying the inequality
constraints as equality constraints.
\item[Ain] Matrix for the inequality constraints \eqn{ {\bf A}_{in}
{\bf p} > {\bf b}_{in}}{A_in p > b}. 
\item[bin] vector in the inequality constraints. 
} 

\end{ldescription}
\end{Arguments}
\begin{Details}\relax
This solves the problem:

\deqn{ minimise~ \| { \bf W}^{1/2} ({ \bf Xp - y} ) \|^2  +  \sum_{i=1}^m
\lambda_i {\bf p^\prime S}_i{\bf p} }{minimise || W^0.5 (Xp-y) ||^2 + 
lambda_1 p'S_1 p + lambda_1 p'S_2 p + . . .}
subject to constraints \eqn{ {\bf Cp}={\bf c}}{Cp=c} and \eqn{ {\bf
A}_{in}{\bf p}>{\bf b}_{in}}{A_in p > b_in}, w.r.t. \eqn{\bf p}{p} given the
smoothing parameters \eqn{\lambda_i}{lambda_i}.
\eqn{ {\bf X}}{X} is a design matrix, \eqn{\bf p}{p} a parameter vector, 
\eqn{\bf y}{y} a data vector, \eqn{\bf W}{W} a diagonal weight matrix,
\eqn{ {\bf S}_i}{S_i} a positive semi-definite matrix  of coefficients
defining the ith penalty and \eqn{\bf C}{C} a matrix of coefficients 
defining the linear equality constraints on the problem. The smoothing
parameters are the \eqn{\lambda_i}{lambda_i}. Note that \eqn{ {\bf X}}{X}
must be of full column rank, at least when projected  into the null space
of any equality constraints. \eqn{ {\bf A}_{in}}{A_in} is a matrix of
coefficients defining the inequality constraints, while \eqn{ {\bf
b}_{in}}{b_in} is a vector involved in defining the inequality constraints.  

Quadratic programming is used to perform the solution. The method used
is designed for maximum stability with least squares problems:
i.e. \eqn{ {\bf X}^\prime {\bf X}}{X'X} is not formed explicitly. See
Gill et al. 1981.
\end{Details}
\begin{Value}
The function returns an array containing the estimated parameter
vector.
\end{Value}
\begin{Author}\relax
Simon N. Wood \email{simon@stats.gla.ac.uk}
\end{Author}
\begin{References}\relax
Gill, P.E., Murray, W. and Wright, M.H. (1981) Practical Optimization. Academic
Press, London. 

Wood, S.N. (1994) Monotonic smoothing splines fitted by cross validation SIAM
Journal on Scientific Computing 15(5):1126-1133

\url{http://www.stats.gla.ac.uk/~simon/}
\end{References}
\begin{SeeAlso}\relax
\code{\LinkA{mgcv}{mgcv}} \code{\LinkA{mono.con}{mono.con}}
\end{SeeAlso}
\begin{Examples}
\begin{ExampleCode}
# first an un-penalized example - fit E(y)=a+bx subject to a>0
set.seed(0)
n<-100
x<-runif(n);y<-x-0.2+rnorm(n)*0.1
M<-list(X=matrix(0,n,2),p=c(0.1,0.5),off=array(0,0),S=list(),
Ain=matrix(0,1,2),bin=0,C=matrix(0,0,0),sp=array(0,0),y=y,w=y*0+1)
M$X[,1]<-1;M$X[,2]<-x;M$Ain[1,]<-c(1,0)
pcls(M)->M$p
plot(x,y);abline(M$p,col=2);abline(coef(lm(y~x)),col=3)

# Penalized example: monotonic penalized regression spline .....

# Generate data from a monotonic truth.
x<-runif(100)*4-1;x<-sort(x);
f<-exp(4*x)/(1+exp(4*x));y<-f+rnorm(100)*0.1;plot(x,y)
dat<-data.frame(x=x,y=y)
# Show regular spline fit (and save fitted object)
f.ug<-gam(y~s(x,k=10,bs="cr"));lines(x,fitted(f.ug))
# Create Design matrix, constraints etc. for monotonic spline....
sm<-smooth.construct(s(x,k=10,bs="cr"),dat,knots=NULL)
F<-mono.con(sm$xp);   # get constraints
G<-list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,p=sm$xp,y=y,w=y*0+1)
G$Ain<-F$A;G$bin<-F$b;G$S<-sm$S;G$off<-0

p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

fv<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv,col=2)

# now a tprs example of the same thing....

f.ug<-gam(y~s(x,k=10));lines(x,fitted(f.ug))
# Create Design matrix, constriants etc. for monotonic spline....
sm<-smooth.construct(s(x,k=10,bs="tp"),dat,knots=NULL)       
xc<-0:39/39 # points on [0,1]  
nc<-length(xc)  # number of constraints
xc<-xc*4-1  # points at which to impose constraints
A0<-Predict.matrix(sm,data.frame(x=xc)) 
# ... A0
A1<-Predict.matrix(sm,data.frame(x=xc+1e-6)) 
A<-(A1-A0)/1e-6    
# ... approx. constraint matrix (A%*%p is -ve spline gradient at points xc)
G<-list(X=sm$X,C=matrix(0,0,0),sp=f.ug$sp,y=y,w=y*0+1,S=sm$S,off=0)
G$Ain<-A;    # constraint matrix
G$bin<-rep(0,nc);  # constraint vector
G$p<-rep(0,10);G$p[10]<-0.1  
# ... monotonic start params, got by setting coefs of polynomial part
p<-pcls(G);  # fit spline (using s.p. from unconstrained fit)

fv2<-Predict.matrix(sm,data.frame(x=x))%*%p
lines(x,fv2,col=3)
\end{ExampleCode}
\end{Examples}

