qle: Simulated quasi-likelihood parameter estimation

In qle: Simulation-Based Quasi-Likelihood Estimation

Description

This is the main estimation function of the simulated quasi-likelihood estimation approach.

Usage

qle(qsd, sim, ..., nsim, x0 = NULL, obs = NULL, Sigma = NULL,
  global.opts = list(), local.opts = list(), method = c("qscoring",
  "bobyqa", "direct"), qscore.opts = list(), control = list(),
  errType = "kv", pl = 0, use.cluster = FALSE, cl = NULL,
  iseed = NULL, plot = FALSE)

Arguments

qsd	object of class QLmodel
sim	user-defined simulation function, see details
...	further arguments passed to 'sim'
nsim	optional, number of simulation replications at each new sample point, set by 'qsd$nsim' (default)
x0	optional, numeric vector of starting parameters
obs	optional, numeric vector of observed statistics, overwrites 'qsd$obs' if given
Sigma	optional, constant variance matrix estimate of statistics (see details)
global.opts	options for global search phase
local.opts	options for local search phase
method	vector of names of local search methods which are applied in consecutive order
qscore.opts	list of control arguments passed to qscoring
control	list of control arguments passed to any of the routines defined in 'method'
errType	type of prediction variances, choose one of "kv","cv", "max" (see details)
pl	print level, use pl>0 to print intermediate results
use.cluster	logical, FALSE (default), whether to use the cluster environment 'cl' for computations other than model simulations or a multicore forking which requires to set options(qle.multicore="mclapply") using at least a number of cpus cores >1, e.g. options(mc.cores=2L).
cl	cluster object, NULL (default), of class "MPIcluster", "SOCKcluster", "cluster"
iseed	integer, seed number, NULL (default) for default seeding of the random number generator (RNG) stream for each worker in the cluster or for parallel processing by "mclapply", if available on non windows platforms. Note that only using the RNG L'Ecuyer-CMRG leads to reproducible results. Only for iseed different from NULL a seed is set including any cluster worker.
plot	if TRUE, plot newly sampled points (for 2D-parameter estimation problems only)

Details

The function sequentially estimates the unknown model parameter. Basically, the user supplies a simulation function 'sim' which must return a vector of summary statistics (as the outcome of model simulations) and expects a vector of parameters as its first input argument. Further arguments can be passed to the simulation function by the '...' argument. The object 'qsd' aggregates the type of variance matrix approximation, the data frame of simulation runs, the initial sample points and the covariance models of the involved statistics (see QLmodel). In addition, it sets the criterion function by 'qsd$criterion', which is either used to monitor the sampling process or minimized itself. The user also has the option to choose among different types of prediction variances: either "kv" (kriging variances), "cv" (cross-validation-based variances) or the maximum of both, set by "max", are available.

Criterion functions

The QD as a criterion function follows the quasi-likelihood estimation principle (see vignette) and seeks a solution to the quasi-score equation. Besides, the Mahalanobis distance (MD) as an alternative criterion function has a more direct interpretation. It can be seen as a (weighted or generalized) least squares criterion depending on the employed type of variance matrix approximation. For this reason, we support several types of variance matrix approximations. In particular, given 'Sigma' and setting 'qsd$var.type' equal to "const" treats 'Sigma' as a constant estimate throughout the whole estimation procedure. Secondly, if 'Sigma' is supplied and used as an average variance approximation (see covarTx), it is considered an initial variance matrix approximation and recomputed each time an approximate (local) minimizer of the criterion function is found. This is commonly known as an iterative update strategy of the variance matrix. Opposed to this, setting 'qsd$var.type' equal to "kriging" corresponds to continuously updating the variance matrix each time a new criterion function value is required at any point of the parameter space. In this way the algorithm can also be seen as a simulated version of a least squares approach or even as a special case of the simulated method of moments (see, e.g. [3]). Note that some input combinations concerning the variance approximation types are not applicable since the criterion "qle", which uses the QD criterion function, is not applicable to a constant variance matrix approximation.

Monte Carlo (MC) hypothesis testing

The algorithm sequentially evaluates promising local minimizers of the criterion function during the local phase in order to assess the plausibility of being an approximate root of the corresponding quasi-score vector. We use essentially the same MC test procedure as in qleTest. First, having found a local minimum of the test statistic, i.e. the criterion function, given the data, new observations are simulated w.r.t. to the local minimizer and the algorithm re-estimates the approximate roots for each observation independently. If the current minimizer is accepted as an approximate root at some significance level 'local.opts$alpha', then the algorithm stays in its local phase and continues sampling around the current minimizer accoring to its asymptotic variance (measured by the inverse of the predicted quasi-information) and uses the additional simulations to improve the current kriging approximations. Otherwise we switch to the global phase and do not consider the current minimizer as an approximate root.

This procedure also allows for a stopping condition derived from the reults of the MC test. We can compare the estimated mean squared error (MSE) with the predicted error of the approximate root by its relative difference and terminate in case this value drops below a user-defined bound 'perr_tol' (either as a scalar value or numeric vector of length equal to the dimension of the unknown parameter). A value close to zero suggests a good match of both error measures. The testing procedure is disabled by default. Set 'local.opts$test=TRUE' for testing an approximate root during the estimation. In this case, if for some value of the criterion function its value exceeds 'local.opts$ftol_abs', then the corresponding minimizer is tested for an approximate root. Otherwise the last evaluation point is used as a starting point for next local searches by a multistart approach when the algorithm is in its global phase. Note that this approach has the potential to escape regions where the criterion function value is quite low but, however, is not considered trustworthy as an approximate root. If testing is disabled, then a local minimizer whose criterion function dropds below the upper bound 'local.opts$ftol_abs' is considered as a root of the quasi-score and the algorithm stays in or switches to its local phase. The same holds, if the quasi-score vector matches the numerical tolerance for being zero. A multistart approach for searching a minimizer during the local phase is only enabled if any of the local search methods fail to converge since most of the time the iteration is started from a point which is already near a root. The re-estimation of parameters during the test procedure also uses this type of multistart approach and cannot be changed by the user.

If one of the other termination criteria is met in conjunction with a neglectable value of the criterion function, we say that the algorithm successfully terminated and converged to a local minimizer of the criterion function which could be an approximate root of the quasi-score vector. This can be further investigated by a goodness-of-fit test in order to assess its plausibility (see qleTest) and quantify the empirical and predicted estimation error of the parameter. If we wish to improve the final estimate the algorithm allows for a simple warm start strategy though not yet as an fully automated procedure. A restart can be based on the final result of the preceeding run. We only need to extract the object 'OPT$qsd' and pass it as an input argument to function qle.

Sampling new points

The QL estimation approach dynamically switches from a local to a global search phase and vise versa for sampling new promising candidates for evaluation, that is, performing new simulations of the statistical model. Depending on the current value of the criterion function three different sampling criteria are used to select next evaluation points which aim on potentially improving the quasi-score or criterion function approximation. If a local minimizer of the criterion function has been accepted as an approximate root, then a local search tries to improve its accuracy. The next evaluation point is either selected according to a weighted minimum-distance criterion (see [2] and vignette), for the choice 'nextSample' equal to "score", or by maximizing the weighted variance of the quasi-score vector in case 'nextSample' is equal to "var". In all other cases, for example, if identifiable roots of the QS could not be found or the (numerical) convergence of the local solvers failed, the global phase of the algorithm is invoked and selects new potential candidates accross the whole search space based on a weighted selection criterion. This assigns large weights to candidates with low criterion function values and vise versa. During both phases the cycling between local and global candidates is controlled by the weights 'global.opts$weights' and 'locals.opts$weights', respectively. Besides this, the smaller the weights the more the candidates tend to be globally selected and vice versa during the global phase. Within the local phase, weights approaching one result in selecting candidates close to the current minimizer of the criterion function. Weights approximately zero maximize the minimum distance between candidates and previously sampled points and thus densely sample the search space almost everywhere if the algorithm is allowed to run infinitely. The choice of weights is somewhat ad hoc but may reflect the users preference on guiding the whole estimation more biased towards either a local or global search. In addition the local weights can be dynamically adjusted if 'useWeights' is FALSE depending on the current progress of estimation. In this case the first weight given by 'locals.opts$weights' is initially used for this kind of adjustment. Make sure to export all functions to the cluster environment 'cl' beforehand, loading required packages on each cluster/compute node, which are used in the model simulation function (see clusterExport and clusterApply).

Parallel options

Parallel processing of all computations is supported either by the multicore approach (spawning/forking tasks by "mclapply" for non Windows-based systems) using the parallel package or a (user-defined) cluster object, e.g. created by makeCluster, currently of class "MPIcluster","SOCKcluster","cluster" which supports calling the function parLapply. By default there is no parallel processing. A cluster object will be automatically created for functions which also could take such object as an argument (if NULL) locally on a single host and finalized after function termination. In this case, the cluster is set up based on a local forking (under Linux) or as a local PSOCK connection with a number of CPUs available for other operating systems. The option "mc.cores", e.g. options(mc.cores=2L, sets the number of cores used for both types of parallel computations. One can also set an integer value 'iseed' as a seed to initialize each worker/thread, see also clusterSetRNGStream, for reproducible results of any function involving random outputs if mc.cores>1 and RNG kind L'Ecuyer-CMRG. Seeding is either done via setting iseed once the user calls the function qle or beforehand using a cluster defined elsewhere in the user calling program. Then the user should set iseed=NULL in order to not reinitialize the seed. The seed passed is stored in the attribute 'optInfo$iseed' of the final return value. If no cluster object is provided, the user sets mc.cores>1L and options(qle.multicore="mclapply"), then most computations and model simulations are performed by function "mclapply" on a single host. In particular, the user-defined stochastic model simulations can be performed in a cluster environment 'cl' whereas all other computations are spawned to the cpus of a single host setting use.cluster=FALSE. We recommend to specify a cluster object once in the user calling program and pass it to functions which take an cluster object as their argument, e.g. qle. Auxiliary computations, e.g. local optimizations or root findings by multiSearch, should be run in parallel on a single host with many CPUs, e.g. setting options(qle.multicore="mclapply") and mc.cores>1. For parallel computations without using a cluster object on a local host simply set both options to "mclapply" and mc.cores>1.

For a 2D parameter estimation problem the function can visualize the sampling and selection process, which requires an active 2D graphical device in advance.

The following controls 'local.opts' for the local search phase are available:

• ftol_rel: upper bound on relative change in criterion function values

• lam_max: upper bound on the maximum eigenvalue of the generalized eigenvalue decomposition of the quasi-information matrix and estimated interpolation error (variance) of quasi-score. This stops the main iteration sampling new locations following the idea that in this case the quasi-score interpolation error has dropped below the estimated precision at best measured by quasi-information matrix for 'global.opts$NmaxLam' consecutive iterates.

• pmin: minimum required probability that a new random candidate point is inside the parameter space. Dropping below this value triggers the global phase. This might indicate that the inverse of the quasi-information matrix does not sufficiently reflect the variance of the current parameter estimate due to a sparse sample or the (hyper)box constraints of the parameter space could be too restrictive.

• nsample: sampling size of candidate locations at the local phase

• weights: vector of weights, 0≤q\code{weights}≤q 1, for local sampling

• useWeights: logical, if FALSE (default), dynamically adjust the weights, see vignette

• ftol_abs: upper bound on the function criterion: values smaller trigger the local phase treating the current minimzer as an approximate root otherwise forces the algorithm to switch to the global phase and vice versa unless testing is enabled in which case, depending on the result of testing, could even stay in the local phase.

• eta: values for decrease and increase of the local weights for the sampling criterion, which is intended to faciliate convergence while sampling new points more and more around the current best parameter estimate.

• alpha: significance level for computation of empirical quantiles of one of the test statistics, that is, testing a parameter to be a root of the quasi-score vector in probability.

• perr_tol upper bound on the relative difference of the empirical and predicted error of an approximate root

• nfail: maximum number of consecutive failed iterations

• nsucc: maximum number of consecutive successful iterations

• nextSample: either "score" (default) or "var" (see details)

The following controls 'global.opts' for the global search phase are available:

• stopval: stopping value related to the criterion function value, the main iteration terminates as soon as the criterion function value drops below this value. This might be preferable to a time consuming sampling procedure if one whishes to simply minimize the criterion function or find a first approximation to the unknown model parameter.

• lam_rel: upper bound on the relative change of 'lam_max', see 'local.opts' The algorithm terminates its value drops below after a number of 'global.opts$NmaxCV' consecutive iterations.

• xtol_rel: relative change of found minimizer of the criterion function or root of quasi-score.

• xscale scaling factors for 'x0', default set to rep(1,length(x0)), i.e. typical magnitude of parameter components

• maxiter: maximum allowed global phase iterations

• maxeval: maximum allowed global and local iterations

• sampleTol: minimum allowed distance between sampled locations at global phase

• weights: vector of \code{weights}>0 for global sampling

• nsample: sampling size of candidate locations at the global phase

• NmaxRel: maximum number of consecutive iterates until stopping according to 'xtol_rel'

• NmaxCV: maximum number of consecutive iterates until stopping according to 'lam_rel'

• NmaxSample: maximum number of consecutive iterations until stopping according to 'sampleTol'

• NmaxLam: maximum number of consecutive iterations until stopping for which the generalized eigenvalue of the variance of the quasi-score vector within the kriging approximation model and its total variance measured by the quasi-information matrix at some estimated parameter drops below the upper bound 'local.opts$lam_max'

• NmaxQI: maximum number of consecutive iterations until stopping for which the relative difference of the empirical error and predicted error of an approximate root drops below 'perr_tol'

• Nmaxftol: maximum number of consecutive iterations until stopping for which the relative change in the values of the criterion function drops below 'local.opts$ftol_rel'

• nstart: number of random starting points, 10*(xdim+1) (default), if local searches do not converge

Value

List of the following objects:

par	final parameter estimate
value	value of criterion function
ctls	a data frame with values of the stopping conditions
qsd	final QLmodel object, including all sample points and covariance models
cvm	CV fitted covariance models
why	names of matched stopping conditions
final	final local minimization results of the criterion function, see searchMinimizer
score	quasi-score vector or the gradient of the Mahalanobis distance
convergence	logical, whether the QL estimation has converged, see details

Attributes:

tracklist	an object (list) of class QDtrack containing local minimization results, evaluated sample points and the status of the corresponding iterations
optInfo	a list of arguments related to the estimation procedure:

• x0: starting parameter vector

• W: final weight matrix, equal to quasi-information matrix at theta, used for both the variance average approximation, if applicable, and as the predicted variance for the (local) sampling of new candidate points according to a multivariate normal distribution with this variance and the current root as its mean parameter.

• theta: the parameter corresponding to W, typically an approximate root or local minimzer of the criterion function

• last.global: logical, whether last iteration sampled a point globally

• minimized: whether last local minimization was successful

• useCV: logical, whether the cross-validation approach was applied

• method: name of final search method applied

• nsim: number of simulation replications at each evaluation point

• iseed the seed to initialize the RNG

Author(s)

M. Baaske

See Also

mahalDist, quasiDeviance, qleTest

Examples

data(normal)
 
# main estimation with new evaluations
# (simulations of the statistical model)
OPT <- qle(qsd,qsd$simfn,nsim=10,
		    global.opts=list("maxeval"=1),
 		local.opts=list("test"=FALSE)) 

qle documentation built on May 2, 2019, 9:55 a.m.