| 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 |
w |
weight vector of the same length as ymyhat or NULL |
jacobian |
function which returns the Jacobian matrix of ymyhat,
namely |
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")
printm(betaLS3,betaLS4)
cat("\n\nComparizon of the Results: gradient\n")
printm(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")
printm(betaLS,betaLS2)
cat("\n\nComparizon of the Results: gradient\n")
printm(gradLS,gradLS2)
## End(Not run)