| JacobianMat {lazy.mat} | R Documentation |
This function returns
when f is a scalar valued function,
Jacobian[j] = d f(x) / d x[j] .
otherwise,
Jacobian[i,j] = d f(x)[i] / d x[j]
JacobianMat(..x.., ..func.., ..eps.. = 1e-06, ...)
..x.. |
A vector of parameters |
..func.. |
function which has x as the first argument |
..eps.. |
small value |
... |
additional parameters to "..func.." |
This program calculates the Jacobian matrix of "..func.." with respect to "..x.." where "..func.." is the scalar or vector valued function and "..x.." is a scalar or vector valued parameters.
The additional paramters to "..func.." can be passed
through ... argument.
In addition, the values of those variables not defined in "..func.."
will be obtained from the environment in which "..func.." is defined.
Typical usage of this function is to calculate the Jacobian matrix
of a statistical model where y is approximated by yhat(param)
by the method of least squares where y and yhat are nobs x 1 vectors.
In this case, ..f..(..x..) is the function to calculate yhat(param).
When yhat(param) is multivariate, denoted by nobs x nvar Yhat(param) matrix,
..f..(..x..) must return c(Yhat(param)).
See the example section of NR
and GN functions of this package,
or the examples in here ExamplesOfGNNR.
When ..func.. is a scalar valued function,
Jacobian[j] = d func(x) / d x[j] .
Otherwise
Jacobian[i,j] = d func(x)[i] / d x[j] .
# scalar valued target function w/o additional arguments
# value of a and d are from parent env.
f <- function(x){ a*log(x)^2+d*x+1 }
dfdx <- function(x, a){
cat("\nThe value to be used: a = ", a, ", d = ", d, "\n")
2*a*log(x)/x+d
}
a <- 1; d <- 2; x <- 1
Jac <-JacobianMat( x, f )
Print( Jac, dfdx( x, a ) )
# scalar valued target function with one additional argument
# a is an argument: d is from parent emv
ff <- function(x, a){ a*log(x)^2+d*x+1 }
dffdx <- function(x, a){
cat("\nThe value to be used: a = ", a, ", d = ", d, "\n")
2*a*log(x)/x+d
}
a <- NA; d <- 2; x <- 1
Jac <-JacobianMat( x, ff, a=1 )
Print( Jac, dffdx( x, a=1 ) )
temp <- function( x, a, d ){
# scalar valued target function
ff <- function(x, a){ a*log(x)^2+d*x+1 }
dfdx <- function(x, a){
cat("\nThe value to be used: a = ", a, ", d = ", d, "\n")
2*a*log(x)/x+d
}
# This d will be used when evaluating ff
Print("in temp", d)
Jac <-JacobianMat( x, ff, a=a )
Print( Jac, dfdx( x, a ) )
}
Print("in global", d)
temp( x=1, a=1, d=2 )
temp( x=1, a=1, d=5 )
# vector valued function w/o additional parameters
# a and b are from parent env.
f2 <- function( x ){
f1 <- a[1]*x[1] + a[2]*x[2]^2
f2 <- b[1]*x[1]^2 + b[2]*x[2]
f3 <- a[1]*x[1]^2 + b[1]*x[2]^2
return( c(f1,f2,f3) )
}
# analytic Jacobian
df2dx <- function( x ){
df1dx1 <- a[1] ; df1dx2 <- 2*a[2]*x[2]
df2dx1 <- 2*b[1]*x[1]; df2dx2 <- b[2]
df3dx1 <- 2*a[1]*x[1]; df3dx2 <- 2*b[1]*x[2]
Jac <- matrix(c(df1dx1, df2dx1, df3dx1, df1dx2, df2dx2, df3dx2), 3)
return( Jac )
}
a=1:2; b=2:3; x=1:2
Jac=JacobianMat( x, f2 )
Print(df2dx(x), Jac)
# vector valued function with an additional parameter
ff2 <- function( x, a ){
f1 <- a[1]*x[1] + a[2]*x[2]^2
f2 <- b[1]*x[1]^2 + b[2]*x[2]
f3 <- a[1]*x[1]^2 + b[1]*x[2]^2
return( c(f1,f2,f3) )
}
# analytic Jacobian
dff2dx <- function( x, a ){
df1dx1 <- a[1] ; df1dx2 <- 2*a[2]*x[2]
df2dx1 <- 2*b[1]*x[1]; df2dx2 <- b[2]
df3dx1 <- 2*a[1]*x[1]; df3dx2 <- 2*b[1]*x[2]
Jac <- matrix(c(df1dx1, df2dx1, df3dx1, df1dx2, df2dx2, df3dx2), 3)
return( Jac )
}
a=NA; b=2:3; x=1:2
Jac=JacobianMat( x, ff2, a=1:2 )
Print(dff2dx(x,a=1:2), Jac)
## Not run:
#
# Calculation of numerical derivatives by JacobianMat
#
#
# Choose one from below.
if(0){
# symbolic definition of a Bivariate function of interest 1
symbolic_f <- expression( A*x^3 + B*x*y + C*y^2 + D )
# symbolic differentiation of f
symbolic_df <- deriv(symbolic_f, c("x","y"), hessian=TRUE)
# symbolic definition of a Univariatefunction of interest 1
symbolic_f <- expression( A*x^2 + B*x + D )
# symbolic differentiation of f
symbolic_df <- deriv(symbolic_f, c("x"), hessian=TRUE)
# symbolic definition of a Univariatefunction of interest 2
symbolic_f <- expression( A*sin(x^2)+x )
# symbolic differentiation of f
symbolic_df <- deriv(symbolic_f, c("x"), hessian=TRUE)
} # end of if(0)
# list of constants in f
CONSTANTS_OF_FUNC=list(A=1, B=1, C=1, D=1)
# Convert the symbolic function to numeric one.
#
# numeric version of the function and its first and second derivatives
f <- function(x,y,z){eval(symbolic_f, envir=CONSTANTS_OF_FUNC)}
df <- function(x,y,z){
g=attributes(
do.call( function( ... ){ eval( symbolic_df, envir=CONSTANTS_OF_FUNC ) }
, as.list(x,y,z) )
)$gradient
cn=colnames(g); g=c(g); names(g)=cn
return(g)
}
df2 <- function(x,y,z){
H=attributes(
do.call( function( ... ){ eval( symbolic_df, envir=CONSTANTS_OF_FUNC ) }
, as.list(x,y,z) )
)$hessian
dn=dimnames(H)[-1]; dim(H)=dim(H)[-1]; dimnames(H)=dn
return( H )
}
# vector argument version
ff <- function( argvec ){
return( do.call(f, as.list(argvec)) )
}
dff <- function( argvec ){
return( do.call( df, as.list(argvec) ) )
}
dff2 <- function( argvec ){
return( do.call( df2, as.list(argvec) ) )
}
#
# Comparison of symbolic and numerical differentiation
#
# Choose one from below.
if(0){
#
# Univariate case
#
x=seq(-1,1,0.1)
for( i in 1:length(x) ){
# true value
d1=df( x[i] )
d2=df2( x[i] )
# numerical differentiation
D1=JacobianMat( x[i], f )
D2=JacobianMat( x[i], function(x){JacobianMat( x, f, ..eps..=1e-4 )})
Print(x[i], d1, D1, d2, D2)
# Print(x[i], d1, D1, d1-D1, d2, D2, d2-D2, fmt="7.5")
}
#
# Bivariate case
#
x=seq(-1,1,0.5)
y=seq(-1,1,0.5)
for( i in 1:length(x) ){
for( j in 1:length(x) ){
xy=c(x[i], y[j])
# true value
d1=dff( xy )
d2=dff2( xy )
# numerical differentiation
D1=JacobianMat( xy, ff )
D2=JacobianMat( xy, function(x){JacobianMat( x, ff, ..eps..=1e-4 )})
Print(xy, d1, D1, d2, D2)
# Print(x[i], d1, D1, d1-D1, d2, D2, d2-D2, fmt="7.5")
}
}
} # end of if(0)
## End(Not run) # end of dontrun