Calculation of Jacobian Matrix

Description

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]

Usage

JacobianMat(..x.., ..func.., ..eps.. = 1e-06, ...)

Arguments

Details

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.

Value

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] .

Examples

# 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