GN {lazy.mat}

R Documentation

Gauss-Newton Method for Least Squares Model Fitting

Description

Gauss-Newton Method for Least Squares Model Fitting

Usage

GN(
  x,
  ymyhat,
  w = NULL,
  jacobian = NULL,
  sd = 0,
  lambda = 1e-04,
  maxiter = 100,
  eps = 1e-09,
  epsg = 1e-09,
  epsx = 1e-09,
  normalize = NULL,
  check_pm = 0,
  maxiter2 = 20,
  print = 1,
  ...
)

Arguments

x	initial parameter vector
ymyhat	residual vector of the form y - yhat(x).
w	weight vector of the same length as ymyhat or NULL
jacobian	function which returns the Jacobian matrix of ymyhat, namely - d yhat(x) / d x.
sd	= 1 to use steepest descent
lambda	constant for quasi Levenberg-Marquardt
maxiter	max # of iterations
eps	eps for the relative change of the function value
epsg	eps for the mazimum absolute value of the gradient
epsx	eps for the maximum absolute change of the parameter value
normalize	a function to normalize the solution after each iteration
check_pm	= 1 to verify partial match of the ... argument
maxiter2	max # of stepsize halving
print	= 1 to print the result, = 2 to print history
...	additional parameters to ymyhat and jacobian

Details

parameter values will be updated as

 x_k+1 = x_k - alpha inv(t(w*jacoban) %*% jacobian + lambda I ) %*%
 (w*grad)

where

 grad = sign 2 t(jacoian) %*% w*ymyhat

where sign=1 if jacobian=NULL and sign=-1 if jacobian is given as #' the Jacobian of y-yhat(x),
alpha is the stepsize to be halved, from 1, when obj increases.

When jacobian is NULL, JacobianMat of this package will be used to calculate it numerically.
See the examples in here ExamplesOfGNNR.

Warning
When the formal arguments of this function partially match the names given as ..., this function does NOT work. In order to avoid this, if any of ... matches the formal arguments, specify those formal arguments by their complete names when calling this function.

Value

A list of:
par the parameter value at minimum
objective the function value minimized
grad the gradient vector
hist a matrix containing the iteration history
iterations number of iterations performed

Examples

#
# Fitting logistic function to observed probability.
#

# functions to be used
Prob <- function( beta, X ){
   # probability (logistic function with matrix X, no 1.7)
   if( is.matrix(X) | is.data.frame(X) ) z=X%*%beta
   else z=X*beta
   P=1/(1+exp(-z))
   return( P )
} # end of P


ymyhat <- function( beta, X=NULL, y=NULL ){
   # returns residuals as y-yhat
   resid=y-Prob(beta,X)
   return( resid )
} # end of ymyhat

dyhat <- function( beta, X=NULL, y=NULL ){
   # Jacobian mat:  d yhat_i / d xi_j
   p=Prob(beta,X)
   F=c(p*(1-p))*X
   return( F )
} # end of dyhat


drss <- function( beta, X=NULL, y=NULL ){
   # gradient of rss
   ymyhat=y-Prob(beta,X)
   F=dyhat( beta, X=X, y=y )
   dd=-2*t(F)%*%ymyhat
   return(dd)
} # end of drss


# seed for random numbers
set.seed(1701)

# true value of beta
beta=c(-1,2,3)
# # of parameters
nq=length(beta)

# # of trials per obs
ff=c(1,2,3) # will be recycled.

# generate regressor matrix X
n=50
# regressor matrix as 1 and normal random variable
X=cbind(1, matrix(rnorm(n*(nq-1)),n,nq-1) )
colnames(X)=paste("x",(1:nq)-1,sep="")

# true response probability
p=Prob(beta,X)
# generate f=# of trials, r=# of successes, y=sample proportion
f=sample( ff, n, rep=1 )
r=mapply( function(size,prob) rbinom(1,size,prob) , f, p )
y=r/f



# initial value
beta0=c(3,2,1)
nq=length(beta)

# LS by GN w/0 Jacobian
resGN3=GN( beta0, ymyhat=ymyhat, X=X, y=y )
betaLS3=resGN3$par
rssLS3=resGN3$objective
gradLS3=drss( betaLS3, X=X, y=y )


# LS by GN with Jacobian
resGN4=GN( beta0, ymyhat=ymyhat, jacobian=dyhat, X=X, y=y )
betaLS4=resGN4$par
rssLS4=resGN4$objective
gradLS4=drss( betaLS4, X=X, y=y )



cat("\n\nComparizon of the Results:  parameters\n")
Print(betaLS3,betaLS4)
cat("\n\nComparizon of the Results:  gradient\n")
Print(gradLS3,gradLS4)



## Not run: 

library(nlsr)
yhatmy <- function( beta, X=NULL, y=NULL ){
   # returns residuals as yhat-y to be used by nlfb
   return( -ymyhat( beta, X=X, y=y ) )
} # end of yhatmy


dyhat2 <- function( beta, X=NULL, y=NULL ){
   # Jacobian mat:  d yhat_i / d xi_j as attribute.  to be used by nlfb
   dyhat=0
   attributes(dyhat)$gradient=dyhat( beta, X=X, y=y )
   return( dyhat )
} # end of dyhat2


# by nlfb
resnlsr=nlfb( beta0, yhatmy, X=X, y=y
              , control=list(femax=500, japprox="jacentral") )
betaLS=resnlsr$coeff
gradLS=drss( betaLS, X=X, y=y )


# by nlfb with Jacobian function
resnlsr2=nlfb( beta0, yhatmy, jacfn=dyhat2, X=X, y=y
               , control=list(femax=500) )
betaLS2=resnlsr2$coeff
gradLS2=drss( betaLS2, X=X, y=y )


cat("\n\nComparizon of the Results:  parameters\n")
Print(betaLS,betaLS2)
cat("\n\nComparizon of the Results:  gradient\n")
Print(gradLS,gradLS2)

## End(Not run)

---