[Georob-commits] r14 - in pkg: . R man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 6 14:05:28 CEST 2013
Author: papritz
Date: 2013-09-06 14:05:27 +0200 (Fri, 06 Sep 2013)
New Revision: 14
Modified:
pkg/ChangeLog
pkg/DESCRIPTION
pkg/NAMESPACE
pkg/R/georob.exported.functions.R
pkg/R/georob.private.functions.R
pkg/man/georob.Rd
pkg/man/georob.control.Rd
Log:
exclusive use of nleqslv for solving estimating equations, released as version 0.1-2 to CRAN
M pkg/R/georob.exported.functions.R
M pkg/R/georob.private.functions.R
M pkg/DESCRIPTION
M pkg/ChangeLog
M pkg/man/georob.Rd
M pkg/man/georob.control.Rd
M pkg/NAMESPACE
Modified: pkg/ChangeLog
===================================================================
--- pkg/ChangeLog 2013-07-16 07:17:28 UTC (rev 13)
+++ pkg/ChangeLog 2013-09-06 12:05:27 UTC (rev 14)
@@ -136,3 +136,9 @@
* georob.exported.functions.R (georob, georob.control, bbsolve.control): solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
* georob.private.functions.R (compute.estimating.equations, compute.expanded.estimating.equations, estimating.eqations.xihat, estimate.xihat, georob.fit, gradient.negative.restricted.loglikelihood, negative.restr.loglikelihood, prepare.likelihood.calculations): solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
+
+
+2013-09-06 Andreas Papritz <papritz at env.ethz.ch>
+
+* georob.exported.functions.R (georob, georob.control, bbsolve.control): code for solving estimating equations by BBsolve{BB} commented out (released to CRAN as version 0.1-2)
+* georob.private.functions.R (compute.estimating.equations, compute.expanded.estimating.equations, estimating.eqations.xihat, estimate.xihat, georob.fit, gradient.negative.restricted.loglikelihood, negative.restr.loglikelihood, prepare.likelihood.calculations): code for solving estimating equations by BBsolve{BB} commented out (released to CRAN as version 0.1-2)
Modified: pkg/DESCRIPTION
===================================================================
--- pkg/DESCRIPTION 2013-07-16 07:17:28 UTC (rev 13)
+++ pkg/DESCRIPTION 2013-09-06 12:05:27 UTC (rev 14)
@@ -1,14 +1,14 @@
Package: georob
Type: Package
Title: Robust Geostatistical Analysis of Spatial Data
-Version: 0.1-1
+Version: 0.1-2
Date: 2013-06-20
Authors at R: c(
person( "Andreas", "Papritz", role = c( "cre", "aut" ),
email = "andreas.papritz at env.ethz.ch" ),
person( "Cornelia", "Schwierz", role = "ctb" ))
Depends: R(>= 2.14.0), lmtest, nlme, robustbase, sp(>= 0.9-60)
-Imports: BB, constrainedKriging(>= 0.2-1), nleqslv, quantreg,
+Imports: constrainedKriging(>= 0.2-1), nleqslv, quantreg,
RandomFields(>= 2.0.55), spatialCovariance(>= 0.6-4)
Suggests: geoR
Description: The georob package provides functions for fitting linear models
Modified: pkg/NAMESPACE
===================================================================
--- pkg/NAMESPACE 2013-07-16 07:17:28 UTC (rev 13)
+++ pkg/NAMESPACE 2013-09-06 12:05:27 UTC (rev 14)
@@ -1,6 +1,6 @@
import( stats, parallel )
-importFrom( BB, BBsolve )
+# importFrom( BB, BBsolve )
importFrom( constrainedKriging, covmodel, f.point.block.cov, K, preCKrige )
importFrom( lmtest, waldtest, waldtest.default )
importFrom( nlme, fixef, fixed.effects, ranef, random.effects )
@@ -11,7 +11,7 @@
# exported functions
export(
- bbsolve.control, # ok
+# bbsolve.control, # ok
bwd.transf, # ok
compress, # ok
cv, # ok
Modified: pkg/R/georob.exported.functions.R
===================================================================
--- pkg/R/georob.exported.functions.R 2013-07-16 07:17:28 UTC (rev 13)
+++ pkg/R/georob.exported.functions.R 2013-09-06 12:05:27 UTC (rev 14)
@@ -19,7 +19,7 @@
aniso = c( f1 = 1., f2 = 1., omega = 90., phi = 90., zeta = 0. ),
fit.aniso = c( f1 = FALSE, f2 = FALSE, omega = FALSE, phi = FALSE, zeta = FALSE ),
tuning.psi = 2, initial.param = c( "minimize", "exclude", "no" ),
- root.finding = c( "nleqslv", "bbsolve" ),
+ ## root.finding = c( "nleqslv", "bbsolve" ),
control = georob.control( ... ), verbose = 0,
...
)
@@ -80,6 +80,7 @@
## 2013-06-12 AP substituting [["x"]] for $x in all lists
## 2013-07-02 AP new transformation of rotation angles
## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
+ ## 2013-09-06 AP exclusive use of nleqslv for solving estimating equations
## check whether input is complete
@@ -318,7 +319,7 @@
aniso.missing <- missing( aniso ) && missing( fit.aniso )
initial.param <- match.arg( initial.param )
- root.finding <- match.arg( root.finding )
+ ## root.finding <- match.arg( root.finding )
## compute initial values of variogram and anisotropy parameters
@@ -364,7 +365,7 @@
## estimate model parameters with pruned data set
t.georob <- georob.fit(
- root.finding = root.finding,
+ ## root.finding = root.finding,
slv = TRUE,
envir = envir,
initial.objects = initial.objects,
@@ -396,8 +397,8 @@
zero.dist = control[["zero.dist"]],
nleqslv.method = control[["nleqslv"]][["method"]],
nleqslv.control = control[["nleqslv"]][["control"]],
- bbsolve.method = control[["bbsolve"]][["method"]],
- bbsolve.control = control[["bbsolve"]][["control"]],
+ ## bbsolve.method = control[["bbsolve"]][["method"]],
+ ## bbsolve.control = control[["bbsolve"]][["control"]],
optim.method = control[["optim"]][["method"]],
optim.lower = control[["optim"]][["lower"]],
optim.upper = control[["optim"]][["upper"]],
@@ -434,7 +435,7 @@
## estimate model parameters by minimizing sum( gradient^2)
t.georob <- georob.fit(
- root.finding = root.finding,
+ ## root.finding = root.finding,
slv = FALSE,
envir = envir,
initial.objects = initial.objects,
@@ -466,8 +467,8 @@
zero.dist = control[["zero.dist"]],
nleqslv.method = control[["nleqslv"]][["method"]],
nleqslv.control = control[["nleqslv"]][["control"]],
- bbsolve.method = control[["bbsolve"]][["method"]],
- bbsolve.control = control[["bbsolve"]][["control"]],
+ ## bbsolve.method = control[["bbsolve"]][["method"]],
+ ## bbsolve.control = control[["bbsolve"]][["control"]],
optim.method = control[["optim"]][["method"]],
optim.lower = control[["optim"]][["lower"]],
optim.upper = control[["optim"]][["upper"]],
@@ -508,7 +509,7 @@
if( verbose > 0 ) cat( "computing final parameter estimates ...\n" )
r.georob <- georob.fit(
- root.finding = root.finding,
+ ## root.finding = root.finding,
slv = TRUE,
envir = envir,
initial.objects = initial.objects,
@@ -540,8 +541,8 @@
zero.dist = control[["zero.dist"]],
nleqslv.method = control[["nleqslv"]][["method"]],
nleqslv.control = control[["nleqslv"]][["control"]],
- bbsolve.method = control[["bbsolve"]][["method"]],
- bbsolve.control = control[["bbsolve"]][["control"]],
+ ## bbsolve.method = control[["bbsolve"]][["method"]],
+ ## bbsolve.control = control[["bbsolve"]][["control"]],
optim.method = control[["optim"]][["method"]],
optim.lower = control[["optim"]][["lower"]],
optim.upper = control[["optim"]][["upper"]],
@@ -622,7 +623,7 @@
rq = rq.control(),
lmrob = lmrob.control(),
nleqslv = nleqslv.control(),
- bbsolve = bbsolve.control(),
+ ## bbsolve = bbsolve.control(),
optim = optim.control(),
full.output = TRUE
)
@@ -671,7 +672,9 @@
aux.cov.pred.target = aux.cov.pred.target,
min.condnum = min.condnum,
irf.models = c( "DeWijsian", "fractalB", "genB" ),
- rq = rq, lmrob = lmrob, nleqslv = nleqslv, bbsolve = bbsolve, optim = optim,
+ rq = rq, lmrob = lmrob, nleqslv = nleqslv,
+ ## bbsolve = bbsolve,
+ optim = optim,
full.output = full.output
)
@@ -804,23 +807,23 @@
}
## ======================================================================
-bbsolve.control <-
- function(
- bbsolve.method = c( "2", "3", "1" ),
- bbsolve.control = NULL
- )
-{
-
- ## function sets meaningful defaults for selected arguments of function
- ## BBSolve{BB}
-
- ## 2013-07-12 A. Papritz
-
- list(
- method = as.integer( match.arg( bbsolve.method ) ),
- control = bbsolve.control
- )
-}
+## bbsolve.control <-
+## function(
+## bbsolve.method = c( "2", "3", "1" ),
+## bbsolve.control = NULL
+## )
+## {
+##
+## ## function sets meaningful defaults for selected arguments of function
+## ## BBSolve{BB}
+##
+## ## 2013-07-12 A. Papritz
+##
+## list(
+## method = as.integer( match.arg( bbsolve.method ) ),
+## control = bbsolve.control
+## )
+## }
## ======================================================================
optim.control <-
Modified: pkg/R/georob.private.functions.R
===================================================================
--- pkg/R/georob.private.functions.R 2013-07-16 07:17:28 UTC (rev 13)
+++ pkg/R/georob.private.functions.R 2013-09-06 12:05:27 UTC (rev 14)
@@ -2646,265 +2646,265 @@
## ##############################################################################
-compute.expanded.estimating.equations <-
- function(
- allpar,
- slv,
- envir,
- variogram.model, fixed.param, param.name, aniso.name,
- param.tf, bwd.tf, safe.param,
- lag.vectors,
- XX, min.condnum, rankdef.x, yy, TT,
- psi.function, dpsi.function,
- tuning.psi, tuning.psi.nr,
- irwls.initial, irwls.maxiter, irwls.reltol,
- force.gradient,
- expectations,
- verbose
- )
-{
-
- ## function evaluates the robustified estimating equations of
- ## variogram parameters derived from the Gaussian log-likelihood
-
- ## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
-
- ## select xihat and variogram parameters
-
- xihat <- allpar[ 1:NROW(XX) ]
- adjustable.param <- allpar[ -(1:NROW(XX)) ]
+## compute.expanded.estimating.equations <-
+## function(
+## allpar,
+## slv,
+## envir,
+## variogram.model, fixed.param, param.name, aniso.name,
+## param.tf, bwd.tf, safe.param,
+## lag.vectors,
+## XX, min.condnum, rankdef.x, yy, TT,
+## psi.function, dpsi.function,
+## tuning.psi, tuning.psi.nr,
+## irwls.initial, irwls.maxiter, irwls.reltol,
+## force.gradient,
+## expectations,
+## verbose
+## )
+## {
+##
+## ## function evaluates the robustified estimating equations of
+## ## variogram parameters derived from the Gaussian log-likelihood
+##
+## ## 2013-07-12 AP solving estimating equations by BBsolve{BB} (in addition to nleqlsv)
+##
+## ## select xihat and variogram parameters
+##
+## xihat <- allpar[ 1:NROW(XX) ]
+## adjustable.param <- allpar[ -(1:NROW(XX)) ]
+##
+## ## get lik.item
+##
+## lik.item <- prepare.likelihood.calculations(
+## envir,
+## adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
+## param.tf, bwd.tf, safe.param,
+## lag.vectors,
+## XX, min.condnum, rankdef.x, yy, TT, xihat,
+## psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
+## irwls.initial, irwls.maxiter, irwls.reltol,
+## compute.xihat = FALSE, compute.Q = FALSE,
+## verbose
+## )
+##
+## ## check whether generalized covariance matrix is positive definite
+##
+## if( lik.item[["Valpha"]][["error"]] ) {
+## if( verbose > 0 ) cat(
+## "\n(generalized) correlation matrix Valpha is not positive definite\n"
+## )
+## t.result <- rep( Inf, length( adjustable.param ) )
+## names( t.result ) <- names( adjustable.param )
+## return( t.result )
+## }
+##
+## ## check whether estimating equations should be computed for fixed parameters
+##
+## if( length( adjustable.param ) == 0 && force.gradient ){
+## adjustable.param <- fixed.param
+## }
+##
+## ## evaluate estimating equations
+##
+## ## compute auxiliary items
+##
+## TtT <- as.vector( table( TT ) )
+##
+## ## compute Cov[bhat]
+##
+## r.cov <- compute.covariances(
+## Valpha.objects = lik.item[["Valpha"]],
+## Aalpha = lik.item[["effects"]][["Aalpha"]],
+## Palpha = lik.item[["effects"]][["Palpha"]],
+## rweights = lik.item[["effects"]][["rweights"]],
+## XX = XX, TT = TT, names.yy = names( yy ),
+## nugget = lik.item[["param"]]["nugget"],
+## eta = lik.item[["eta"]],
+## expectations = expectations,
+## cov.bhat = TRUE, full.cov.bhat = TRUE,
+## cov.betahat = FALSE,
+## cov.bhat.betahat = FALSE,
+## cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
+## cov.delta.bhat.betahat = FALSE,
+## cov.ehat = FALSE, full.cov.ehat = FALSE,
+## cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
+## aux.cov.pred.target = FALSE,
+## extended.output = FALSE,
+## verbose = verbose
+## )
+##
+## if( r.cov[["error"]] ) {
+## if( verbose > 0 ) cat(
+## "\nan error occurred when computing the covariances of fixed and random effects\n"
+## )
+## t.result <- rep( Inf, length( adjustable.param ) )
+## names( t.result ) <- names( adjustable.param )
+## return( t.result )
+## }
+##
+## ## estimating equations for xihat
+##
+## eeq.xihat <- estimating.eqations.xihat(
+## res = lik.item[["effects"]][["residuals"]],
+## TT = TT, xihat = xihat,
+## nugget = lik.item[["param"]]["nugget"],
+## eta = lik.item[["eta"]],
+## Valpha.inverse.Palpha = lik.item[["effects"]][["Valpha.inverse.Palpha"]],
+## psi.function = psi.function,
+## tuning.psi = tuning.psi
+## )
+##
+## ## initialize estimating equations for variogram parameters
+##
+## eeq.emp <- rep( NA, length( adjustable.param ) )
+## names( eeq.emp ) <- names( adjustable.param )
+##
+## eeq.exp <- rep( NA, length( adjustable.param ) )
+## names( eeq.exp ) <- names( adjustable.param )
+##
+## ## estimation equation for nugget
+##
+## if( "nugget" %in% names( adjustable.param ) ) {
+##
+## ## compute trace of Cov[ psi( residuals/sqrt(nugget) ) ]
+##
+## eeq.exp["nugget"] <- sum(
+## diag(
+## lik.item[["Valpha"]][["Valpha.inverse"]] %*%
+## ( 1/TtT * lik.item[["Valpha"]][["Valpha.inverse"]] ) %*%
+## r.cov[["cov.bhat"]]
+## )
+## )
+## eeq.emp["nugget"] <- sum(
+## ( lik.item[["effects"]][["z.star"]] )^2 / TtT
+## )
+##
+## }
+##
+## ## estimation equation for spatial nugget
+##
+## if( "snugget" %in% names( adjustable.param ) ) {
+##
+## ## compute trace( Valpha^-1 Cov[bhat] )
+##
+## eeq.exp["snugget"] <- sum(
+## rowSums(
+## (lik.item[["Valpha"]][["Valpha.inverse"]] %*% lik.item[["Valpha"]][["Valpha.inverse"]] ) *
+## r.cov[["cov.bhat"]]
+## )
+## )
+## eeq.emp["snugget"] <- sum( lik.item[["effects"]][["z.star"]]^2 )
+##
+## }
+##
+## ## estimation equation for variance
+##
+## if( "variance" %in% names( adjustable.param ) ) {
+##
+## ## compute trace( Valpha^-1 Cov[bhat] )
+##
+## eeq.exp["variance"] <- sum(
+## rowSums(
+## ( lik.item[["Valpha"]][["Valpha.inverse"]] %*% lik.item[["Valpha"]][["Valpha0"]] %*% lik.item[["Valpha"]][["Valpha.inverse"]] ) *
+## r.cov[["cov.bhat"]]
+## )
+## )
+## eeq.emp["variance"] <- sum(
+## lik.item[["effects"]][["z.star"]] * drop( lik.item[["Valpha"]][["Valpha0"]] %*% lik.item[["effects"]][["z.star"]] )
+## )
+##
+## }
+##
+## ## estimation equations for scale, extra variogram and anisotropy
+## ## parameters
+##
+## extra.par <- names( adjustable.param )[ !(
+## names( adjustable.param ) %in% c( "variance", "snugget", "nugget" )
+## )]
+##
+## for( t.i in extra.par ){
+##
+## ## compute trace( Valpha^-1 * dValpha/dalpha * Valpha^-1 * Cov[bhat] )
+##
+## dValpha <- dcorr.dparam(
+## x = lag.vectors, variogram.model = variogram.model, param = lik.item[["param"]],
+## d.param = t.i,
+## aniso = lik.item[["aniso"]],
+## verbose = verbose
+## )
+## ## if( identical( class( dValpha ), "try-error" ) ){
+## ## if( verbose > 0 ) cat( "error in dcorr.dparam\n\n" )
+## ## t.result <- rep( Inf, length( adjustable.param ) )
+## ## names( t.result ) <- names( adjustable.param )
+## ## return( t.result )
+## ## }
+##
+## eeq.exp[t.i] <- sum(
+## rowSums(
+## (lik.item[["Valpha"]][["Valpha.inverse"]] %*% dValpha %*% lik.item[["Valpha"]][["Valpha.inverse"]]) *
+## r.cov[["cov.bhat"]]
+## )
+## )
+## eeq.emp[t.i] <- sum(
+## lik.item[["effects"]][["z.star"]] * drop( dValpha %*% lik.item[["effects"]][["z.star"]] )
+## )
+##
+## }
+##
+## if( verbose > 1 ) {
+## cat( "\n ",
+## format( c( "min(xihat)", "max(xihat)" ), width = 14, justify = "right" ),
+## "\n", sep =""
+## )
+## cat( " EEQ :",
+## format(
+## signif( range(eeq.xihat), digits = 7 ),
+## scientific = TRUE, width = 14
+## ), "\n", sep = ""
+## )
+## cat( "\n ",
+## format( names( eeq.emp), width = 14, justify = "right" ),
+## "\n", sep =""
+## )
+## cat( " EEQ :",
+## format(
+## signif( eeq.emp / eeq.exp - 1, digits = 7 ),
+## scientific = TRUE, width = 14
+## ), "\n", sep = ""
+## )
+## if( verbose > 2 ){
+## cat( " empirical terms:",
+## format(
+## signif( eeq.emp, digits = 7 ),
+## scientific = TRUE, width = 14
+## ), "\n", sep = ""
+## )
+## cat( " expected terms:",
+## format(
+## signif( eeq.exp, digits = 7 ),
+## scientific = TRUE, width = 14
+## ), "\n", sep = ""
+## )
+## }
+## cat("\n")
+## }
+##
+## ## store terms in lik.item object
+##
+## lik.item[["eeq"]] <- list(
+## eeq.xihat = eeq.xihat,
+## eeq.emp = eeq.emp,
+## eeq.exp = eeq.exp
+## )
+##
+## assign( "lik.item", lik.item, pos = as.environment( envir ) )
+##
+## return( c( eeq.xihat, eeq.emp / eeq.exp - 1. ) )
+##
+## }
- ## get lik.item
-
- lik.item <- prepare.likelihood.calculations(
- envir,
- adjustable.param, variogram.model, fixed.param, param.name, aniso.name,
- param.tf, bwd.tf, safe.param,
- lag.vectors,
- XX, min.condnum, rankdef.x, yy, TT, xihat,
- psi.function, dpsi.function, tuning.psi, tuning.psi.nr,
- irwls.initial, irwls.maxiter, irwls.reltol,
- compute.xihat = FALSE, compute.Q = FALSE,
- verbose
- )
-
- ## check whether generalized covariance matrix is positive definite
-
- if( lik.item[["Valpha"]][["error"]] ) {
- if( verbose > 0 ) cat(
- "\n(generalized) correlation matrix Valpha is not positive definite\n"
- )
- t.result <- rep( Inf, length( adjustable.param ) )
- names( t.result ) <- names( adjustable.param )
- return( t.result )
- }
-
- ## check whether estimating equations should be computed for fixed parameters
-
- if( length( adjustable.param ) == 0 && force.gradient ){
- adjustable.param <- fixed.param
- }
-
- ## evaluate estimating equations
-
- ## compute auxiliary items
-
- TtT <- as.vector( table( TT ) )
-
- ## compute Cov[bhat]
-
- r.cov <- compute.covariances(
- Valpha.objects = lik.item[["Valpha"]],
- Aalpha = lik.item[["effects"]][["Aalpha"]],
- Palpha = lik.item[["effects"]][["Palpha"]],
- rweights = lik.item[["effects"]][["rweights"]],
- XX = XX, TT = TT, names.yy = names( yy ),
- nugget = lik.item[["param"]]["nugget"],
- eta = lik.item[["eta"]],
- expectations = expectations,
- cov.bhat = TRUE, full.cov.bhat = TRUE,
- cov.betahat = FALSE,
- cov.bhat.betahat = FALSE,
- cov.delta.bhat = FALSE, full.cov.delta.bhat = FALSE,
- cov.delta.bhat.betahat = FALSE,
- cov.ehat = FALSE, full.cov.ehat = FALSE,
- cov.ehat.p.bhat = FALSE, full.cov.ehat.p.bhat = FALSE,
- aux.cov.pred.target = FALSE,
- extended.output = FALSE,
- verbose = verbose
- )
-
- if( r.cov[["error"]] ) {
- if( verbose > 0 ) cat(
- "\nan error occurred when computing the covariances of fixed and random effects\n"
- )
- t.result <- rep( Inf, length( adjustable.param ) )
- names( t.result ) <- names( adjustable.param )
- return( t.result )
- }
-
- ## estimating equations for xihat
-
- eeq.xihat <- estimating.eqations.xihat(
- res = lik.item[["effects"]][["residuals"]],
- TT = TT, xihat = xihat,
- nugget = lik.item[["param"]]["nugget"],
- eta = lik.item[["eta"]],
- Valpha.inverse.Palpha = lik.item[["effects"]][["Valpha.inverse.Palpha"]],
- psi.function = psi.function,
- tuning.psi = tuning.psi
- )
-
- ## initialize estimating equations for variogram parameters
-
- eeq.emp <- rep( NA, length( adjustable.param ) )
- names( eeq.emp ) <- names( adjustable.param )
-
- eeq.exp <- rep( NA, length( adjustable.param ) )
- names( eeq.exp ) <- names( adjustable.param )
-
- ## estimation equation for nugget
-
- if( "nugget" %in% names( adjustable.param ) ) {
-
- ## compute trace of Cov[ psi( residuals/sqrt(nugget) ) ]
-
- eeq.exp["nugget"] <- sum(
- diag(
- lik.item[["Valpha"]][["Valpha.inverse"]] %*%
- ( 1/TtT * lik.item[["Valpha"]][["Valpha.inverse"]] ) %*%
- r.cov[["cov.bhat"]]
- )
- )
- eeq.emp["nugget"] <- sum(
- ( lik.item[["effects"]][["z.star"]] )^2 / TtT
- )
-
- }
-
- ## estimation equation for spatial nugget
-
- if( "snugget" %in% names( adjustable.param ) ) {
-
- ## compute trace( Valpha^-1 Cov[bhat] )
-
- eeq.exp["snugget"] <- sum(
- rowSums(
- (lik.item[["Valpha"]][["Valpha.inverse"]] %*% lik.item[["Valpha"]][["Valpha.inverse"]] ) *
- r.cov[["cov.bhat"]]
- )
- )
- eeq.emp["snugget"] <- sum( lik.item[["effects"]][["z.star"]]^2 )
-
- }
-
- ## estimation equation for variance
-
- if( "variance" %in% names( adjustable.param ) ) {
-
- ## compute trace( Valpha^-1 Cov[bhat] )
-
- eeq.exp["variance"] <- sum(
- rowSums(
- ( lik.item[["Valpha"]][["Valpha.inverse"]] %*% lik.item[["Valpha"]][["Valpha0"]] %*% lik.item[["Valpha"]][["Valpha.inverse"]] ) *
- r.cov[["cov.bhat"]]
- )
- )
- eeq.emp["variance"] <- sum(
- lik.item[["effects"]][["z.star"]] * drop( lik.item[["Valpha"]][["Valpha0"]] %*% lik.item[["effects"]][["z.star"]] )
- )
-
- }
-
- ## estimation equations for scale, extra variogram and anisotropy
- ## parameters
-
- extra.par <- names( adjustable.param )[ !(
- names( adjustable.param ) %in% c( "variance", "snugget", "nugget" )
- )]
-
- for( t.i in extra.par ){
-
- ## compute trace( Valpha^-1 * dValpha/dalpha * Valpha^-1 * Cov[bhat] )
-
- dValpha <- dcorr.dparam(
- x = lag.vectors, variogram.model = variogram.model, param = lik.item[["param"]],
- d.param = t.i,
- aniso = lik.item[["aniso"]],
- verbose = verbose
- )
- ## if( identical( class( dValpha ), "try-error" ) ){
- ## if( verbose > 0 ) cat( "error in dcorr.dparam\n\n" )
- ## t.result <- rep( Inf, length( adjustable.param ) )
- ## names( t.result ) <- names( adjustable.param )
- ## return( t.result )
- ## }
-
- eeq.exp[t.i] <- sum(
- rowSums(
- (lik.item[["Valpha"]][["Valpha.inverse"]] %*% dValpha %*% lik.item[["Valpha"]][["Valpha.inverse"]]) *
- r.cov[["cov.bhat"]]
- )
- )
- eeq.emp[t.i] <- sum(
- lik.item[["effects"]][["z.star"]] * drop( dValpha %*% lik.item[["effects"]][["z.star"]] )
- )
-
- }
-
- if( verbose > 1 ) {
- cat( "\n ",
- format( c( "min(xihat)", "max(xihat)" ), width = 14, justify = "right" ),
- "\n", sep =""
- )
- cat( " EEQ :",
- format(
- signif( range(eeq.xihat), digits = 7 ),
- scientific = TRUE, width = 14
- ), "\n", sep = ""
- )
- cat( "\n ",
- format( names( eeq.emp), width = 14, justify = "right" ),
- "\n", sep =""
- )
- cat( " EEQ :",
- format(
- signif( eeq.emp / eeq.exp - 1, digits = 7 ),
- scientific = TRUE, width = 14
- ), "\n", sep = ""
- )
- if( verbose > 2 ){
- cat( " empirical terms:",
- format(
- signif( eeq.emp, digits = 7 ),
- scientific = TRUE, width = 14
- ), "\n", sep = ""
- )
- cat( " expected terms:",
- format(
- signif( eeq.exp, digits = 7 ),
- scientific = TRUE, width = 14
- ), "\n", sep = ""
- )
- }
- cat("\n")
- }
-
- ## store terms in lik.item object
-
- lik.item[["eeq"]] <- list(
- eeq.xihat = eeq.xihat,
- eeq.emp = eeq.emp,
- eeq.exp = eeq.exp
- )
-
- assign( "lik.item", lik.item, pos = as.environment( envir ) )
-
- return( c( eeq.xihat, eeq.emp / eeq.exp - 1. ) )
-
-}
-
## ##############################################################################
negative.restr.loglikelihood <-
@@ -3460,7 +3460,7 @@
georob.fit <-
function(
- root.finding,
+ ## root.finding,
slv,
envir,
initial.objects,
@@ -3489,7 +3489,7 @@
force.gradient,
zero.dist,
nleqslv.method, nleqslv.control,
- bbsolve.method, bbsolve.control,
+ ## bbsolve.method, bbsolve.control,
optim.method, optim.lower, optim.upper, hessian, optim.control,
full.output,
verbose
@@ -3949,7 +3949,7 @@
if( slv ){
- if( identical( root.finding, "nleqslv" ) ){
+ ## if( identical( root.finding, "nleqslv" ) ){
r.root <- nleqslv(
x = c(
@@ -4000,116 +4000,82 @@
r.counts <- c( nfcnt = r.root[["nfcnt"]], njcnt = r.root[["njcnt"]] )
- } else if( identical( root.finding, "bbsolve" ) ) {
+ ## }
+ ##
+ ## else if( identical( root.finding, "bbsolve" ) ) {
+ ##
+ ## r.root <- BBsolve(
+ ## par = c(
+ ## xihat,
+ ## transformed.param[ fit.param ],
+ ## transformed.aniso[ fit.aniso ]
+ ## ),
+ ## fn = compute.expanded.estimating.equations,
+ ## method = bbsolve.method,
+ ## control = bbsolve.control,
+ ## quiet = verbose == 0,
+ ## slv = slv,
+ ## envir = envir,
+ ## variogram.model = variogram.model,
+ ## fixed.param = c(
+ ## transformed.param[ !fit.param ],
+ ## transformed.aniso[ !fit.aniso ]
+ ## ),
+ ## param.name = param.name,
+ ## aniso.name = aniso.name,
+ ## param.tf = param.tf,
+ ## bwd.tf = bwd.tf,
+ ## safe.param = safe.param,
+ ## lag.vectors = lag.vectors,
+ ## XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
+ ## yy = yy, TT = TT,
+ ## psi.function = rho.psi.etc[["psi.function"]],
+ ## dpsi.function = rho.psi.etc[["dpsi.function"]],
+ ## tuning.psi = tuning.psi,
+ ## tuning.psi.nr = tuning.psi.nr,
+ ## irwls.initial = irwls.initial,
+ ## irwls.maxiter = irwls.maxiter,
+ ## irwls.reltol = irwls.reltol,
+ ## force.gradient = force.gradient,
+ ## expectations = expectations,
+ ## verbose = verbose
+ ## )
+ ##
+ ## r.converged <- r.root[["convergence"]] == 0
+ ## r.convergence.code <- r.root[["convergence"]]
+ ## r.counts <- c( nfcnt = r.root[["feval"]], njcnt = NA_integer_ )
+ ##
+ ## r.gradient <- compute.expanded.estimating.equations(
+ ## allpar = r.root[["par"]],
+ ## slv = TRUE,
+ ## envir = envir,
+ ## variogram.model = variogram.model,
+ ## fixed.param = c(
+ ## transformed.param[ !fit.param ],
+ ## transformed.aniso[ !fit.aniso ]
+ ## ),
+ ## param.name = param.name,
+ ## aniso.name = aniso.name,
+ ## param.tf = param.tf,
+ ## bwd.tf = bwd.tf,
+ ## safe.param = safe.param,
+ ## lag.vectors = lag.vectors,
+ ## XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
+ ## yy = yy, TT = TT,
+ ## psi.function = rho.psi.etc[["psi.function"]],
+ ## dpsi.function = rho.psi.etc[["dpsi.function"]],
+ ## tuning.psi = tuning.psi,
+ ## tuning.psi.nr = tuning.psi.nr,
+ ## irwls.initial = irwls.initial,
+ ## irwls.maxiter = irwls.maxiter,
+ ## irwls.reltol = irwls.reltol,
+ ## force.gradient = force.gradient,
+ ## expectations = expectations,
+ ## verbose = verbose
+ ## )
+ ##
+ ## }
- r.root <- BBsolve(
- par = c(
- xihat,
- transformed.param[ fit.param ],
- transformed.aniso[ fit.aniso ]
- ),
- fn = compute.expanded.estimating.equations,
- method = bbsolve.method,
- control = bbsolve.control,
- quiet = verbose == 0,
- slv = slv,
- envir = envir,
- variogram.model = variogram.model,
- fixed.param = c(
- transformed.param[ !fit.param ],
- transformed.aniso[ !fit.aniso ]
- ),
- param.name = param.name,
- aniso.name = aniso.name,
- param.tf = param.tf,
- bwd.tf = bwd.tf,
- safe.param = safe.param,
- lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
- yy = yy, TT = TT,
- psi.function = rho.psi.etc[["psi.function"]],
- dpsi.function = rho.psi.etc[["dpsi.function"]],
- tuning.psi = tuning.psi,
- tuning.psi.nr = tuning.psi.nr,
- irwls.initial = irwls.initial,
- irwls.maxiter = irwls.maxiter,
- irwls.reltol = irwls.reltol,
- force.gradient = force.gradient,
- expectations = expectations,
- verbose = verbose
- )
-#
-# r.root <- BBsolve(
-# par = c(
-# transformed.param[ fit.param ],
-# transformed.aniso[ fit.aniso ]
-# ),
-# fn = compute.estimating.equations,
-# method = bbsolve.method,
-# control = bbsolve.control,
-# quiet = verbose == 0,
-# slv = slv,
-# envir = envir,
-# variogram.model = variogram.model,
-# fixed.param = c(
-# transformed.param[ !fit.param ],
-# transformed.aniso[ !fit.aniso ]
-# ),
-# param.name = param.name,
-# aniso.name = aniso.name,
-# param.tf = param.tf,
-# bwd.tf = bwd.tf,
-# safe.param = safe.param,
-# lag.vectors = lag.vectors,
-# XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
-# yy = yy, TT = TT, xihat = xihat,
-# psi.function = rho.psi.etc[["psi.function"]],
-# dpsi.function = rho.psi.etc[["dpsi.function"]],
-# tuning.psi = tuning.psi,
-# tuning.psi.nr = tuning.psi.nr,
-# irwls.initial = irwls.initial,
-# irwls.maxiter = irwls.maxiter,
-# irwls.reltol = irwls.reltol,
-# force.gradient = force.gradient,
-# expectations = expectations,
-# verbose = verbose
-# )
-#
- r.converged <- r.root[["convergence"]] == 0
- r.convergence.code <- r.root[["convergence"]]
- r.counts <- c( nfcnt = r.root[["feval"]], njcnt = NA_integer_ )
-
- r.gradient <- compute.expanded.estimating.equations(
- allpar = r.root[["par"]],
- slv = TRUE,
- envir = envir,
- variogram.model = variogram.model,
- fixed.param = c(
- transformed.param[ !fit.param ],
- transformed.aniso[ !fit.aniso ]
- ),
- param.name = param.name,
- aniso.name = aniso.name,
- param.tf = param.tf,
- bwd.tf = bwd.tf,
- safe.param = safe.param,
- lag.vectors = lag.vectors,
- XX = XX, min.condnum = min.condnum, rankdef.x = rankdef.x,
- yy = yy, TT = TT,
- psi.function = rho.psi.etc[["psi.function"]],
[TRUNCATED]
To get the complete diff run:
svnlook diff /svnroot/georob -r 14
More information about the Georob-commits
mailing list