[R] How to do a backward calculation for each record in a dataset

[Rolf Turner]

On 02/19/2013 09:44 AM, Berend Hasselman wrote:
> Rolf,
>
> Your attachments got through.
> But where is the function newt(…)?

Ahhhh --- dang!  It's another item from my "personal miscellany package"
which I'm so used to having around that I forgot that other people don't 
have it.
Attached.   Along with its documentation file.

Thanks for pointing out my sin of omission. :-)

     cheers,

         Rolf
-------------- next part --------------
newt <- function(fn,start,...,eps.p = 1e-08,eps.v = NULL,
                 maxit = 50,verb = FALSE)
{
p.o <- start
itno <- 1
repeat {
	fj <- fn(p.o,...)
	v <- fj$fval
	t1 <- if(is.null(eps.v)) NULL else sum(abs(v))
	J <- as.matrix(fj$jacobian)
	if(qr(J)$rank < ncol(J)) {
		cat("Singular Jacobian.\n")
		rslt <- if(is.null(eps.v)) NA else if(t1 < eps.v) p.o
		else NA
		break
	}
	else {
		p.n <- p.o - solve(J) %*% v
		t2 <- max(abs(p.n - p.o))
		if(verb) {
			tmp <- format(round(c(p.o,p.n,v,t2,t1),6))
			np <- length(v)
			v1 <- paste(tmp[1:np],collapse = "  ")
			v2 <- paste(tmp[(np + 1):(2 * np)],collapse = "  ")
			v3 <- paste(tmp[(2 * np + 1):(3 * np)],collapse = "  ")
			v4 <- tmp[3 * np + 1]
			v5 <- tmp[3 * np + 2]
			cat("\nIteration  : ",itno,"\n",sep = "")
			cat("Old par    : ",v1,"\n",sep = "")
			cat("New par    : ",v2,"\n",sep = "")
			cat("Test ch.par: ",v4,"\n",sep = "")
			cat("Fn. vals.  : ",v3,"\n",sep = "")
			if(!is.null(t1))
			  cat("Test f.val: ",v5,"\n",sep = "")
		}
		if((!is.null(t1) && t1 < eps.v) | t2 < eps.p) {
			rslt <- p.n
			break
		}
		itno <- itno + 1
		if(itno > maxit) {
			cat("Newton's method failed to converge in\n")
			cat(maxit,"iterations.\n")
			rslt <- NA
			break
		}
		p.o <- p.n
	}
}
as.vector(rslt)
}
-------------- next part --------------
\name{newt}
\alias{newt}
\title{
Newton's method.
}
\description{
A rather naive general implementation of Newton's method;
but it seems to work.  A lot of the time! [ :-) ]
}
\usage{
newt(fn, start, \dots, eps.p=1e-08, eps.v=NULL,
     maxit=50, verb=FALSE)
}
\arguments{
\item{fn}{
A user-supplied function providing the essential information about
the system of equations to be solved.  The system is assumed to be of
the form \eqn{f(x) = 0} where \eqn{f(x)} is a \eqn{k}-dimensional
function (with components \eqn{f_1(x) \ldots f_k(x)}{f_1(x) ...
f_k(x)} of \eqn{k} variables.

The function \code{fn} must be coded to return a list with components
\code{fval} and \code{jacobian}.  The component \code{fval} must be a
\eqn{k}-dimensional vector (whose \eqn{i^{th}}{i-th} component is the
value of \eqn{f_i(x)}.  The component "jacobian" is the Jacobian of
the function \eqn{f}, i.e. it is a matrix whose
\eqn{(i,j)^{th}}{(i,j)-th} entry is the derivative of \eqn{f_i(x)}
with respect to \eqn{x_j}.
}
\item{start}{
A \eqn{k}-dimensional vector of starting values from which to
iterate toward the solution.
}
\item{...}{
Any auxilliary arguments needed by the function \code{fn}.
}
\item{eps.p}{
The iteration stops if the maximum absolute value of the
change in the parameters \eqn{x_1, \ldots, x_k}{x_1, ..., x_k}
is less than \code{eps.p}.
}
\item{eps.v}{
If this argument is provided the iteration stops if the sum of the
absolute values of the function values \eqn{f_j(x)} is less than
\code{eps.v}.  (Note: If \code{eps.v} is provided then iteration will
cease if EITHER the "\code{eps.p} criterion" or the "\code{eps.v}
criterion" is met.)
}
\item{maxit}{
The maximum number of iterations to attempt before giving
up in disgust.
}
\item{verb}{
Logical scalar; if TRUE a description of the current state
of play is printed out at every iteration.
}}
\value{
If the iterative procedure has converged, to within the specified
tolerance(s), the final value of the k-dimensional vector of
parameter ("x") values.  Otherwise, NA.
}
\details{
If the Jacobian becomes (numerically) singular (as determined by the
qv() function) then the function \code{newt} exits.  If \code{eps.v}
is not provided a value of NA is returned.  If \code{eps.v} IS
provided, and if by some miracle the sum of the absolute values of
the function values \eqn{f_j(x)} is less than \code{eps.v}, then the
current value of \eqn{x} is returned (since this \eqn{x} does satisfy
the set of equations to the specified tolerance).
}
\author{Rolf Turner
  \email{r.turner at auckland.ac.nz}
  \url{http://www.math.unb.ca/~rolf}
}
\examples{
foo <- function(x) {
   fval <- c(x[1]**2 + x[2]**2 - 1, x[2] - x[1])
   jacobian <- matrix(c(2*x[1],2*x[2], -1, 1),byrow=TRUE,ncol=2)
   list(fval=fval,jacobian=jacobian)
}
newt(foo,c(0.5,0.5))   # gives (1/sqrt(2),1/sqrt(2))
newt(foo,c(-0.5,-0.5)) # gives (-1/sqrt(2),-1/sqrt(2))
newt(foo,c(-0.5,0.5))  # Singular Jacobian; NA returned.
newt(foo,c(-0.7,0.5))  # gives (-1/sqrt(2),-1/sqrt(2))
}
\keyword{math}