Smooth.LS {CollocInfer}

Model-Based Smoothing Functions

Description

Perform the inner optimization to estimate coefficients given parameters.

Usage

Smooth.LS(fn,data,times,pars,coefs=NULL,basisvals=NULL,lambda,fd.obj=NULL,
        more=NULL,weights=NULL,quadrature=NULL,likfn = make.id(), 
        likmore = NULL,in.meth='nlminb',control.in,eps=1e-6,
        posproc=FALSE,poslik=FALSE,discrete=FALSE,names=NULL,
        sparse=FALSE)
        
Smooth.multinorm(fn,data,times,pars,coefs=NULL,basisvals=NULL,var=c(1,0.01),
        fd.obj=NULL,more=NULL,quadrature=NULL,in.meth='nlminb',
        control.in,eps=1e-6,posproc=FALSE,poslik=FALSE,discrete=FALSE,
        names=NULL,sparse=FALSE)

Arguments

fn

A function giving the right hand side of a differential/difference equation. The function should have arguments

• times The times at which the RHS is being evaluated.

• x The state values at those times.

• p Parameters to be entered in the system.

• more An object containing additional inputs to fn

It should return a matrix of the same dimension of x giving the right hand side values.

If fn is given as a single function, its derivatives are estimated by finite-differencing with stepsize eps. Alternatively, a list can be supplied with elements:

• fn Function to calculate the right hand side should accept a matrix of state values at .

• dfdx Function to calculate the derivative with respect to x

• dfdp Function to calculate the derivative with respect to p

• d2fdx2 Function to calculate the second derivative with respect to x

• d2fdxdp Function to calculate the second derivative with respect to x and p

These functions take the same arguments as fn and should output multidimensional arrays with the dimensions ordered according to time, state, deriv1, deriv2; here derivatives with respect to x always precede derivatives with respect to p.

data

Matrix of observed data values.

times

Vector observation times for the data.

pars

Initial values of parameters to be estimated processes.

coefs

Vector giving the current estimate of the coefficients in the spline.

basisvals

Values of the collocation basis to be used. This can either be a basis object from the fda package, or a list elements:

• bvals.obs A matrix giving the values of the basis at the observation times

• bvals A matrix giving the values of the basis at the quadrature times

• dbvals A matrix giving the derivative of the basis at the quadrature times

lambda

(Smooth.LS only) Penalty value trading off fidelity to data with fidelity to differential equations.

var

(Smooth.multinorm) A vector of length 2, giving

fd.obj

(Optional) A functional data object; if this is non-null, coefs and basisvals is extracted from here.

more

An object specifying additional arguments to fn.

weights

(Smooth.LS only)

quadrature

Quadrature points, should contain two elements (if not NULL)

• qpts Quadrature points; defaults to midpoints between knots

• qwts Quadrature weights; defaults to normalizing by the length of qpts.

in.meth

Inner optimization function to be used, currently one of 'nlminb', 'MaxNR', 'optim' or 'SplineEst'. The last calls SplineEst.NewtRaph. This is fast but has poor convergence.

control.in

Control object for inner optimization function.

eps

Finite differencing step size, if needed.

posproc

Should the state vector be constrained to be positive? If this is the case, the state is represented by an exponentiated basis expansion in the proc object.

poslik

Should the state be exponentiated before being compared to the data? When the state is represented on the log scale (posproc=TRUE), this is an alternative to taking the log of the data.

discrete

Is this a discrete or continuous-time system?

names

The names of the state variables if not given by the column names of coefs.

sparse

Should sparse matrices be used for basis values? This option can save memory when ProfileGausNewt and SplineEstNewtRaph are called. Otherwise sparse matrices will be converted to full matrices and this can slow the code down.

likfn

Defines a map from the trajectory to the observations. This should be in the same form as fn. If a function is given, derivatives are estimated by finite differencing, otherwise a list is expected to provide the same derivatives as fn. If poslik=TRUE, the states are exponentiated before the likfn is evaluated and the derivatives are updated to account for this. Defaults to the identity transform.

likmore

A list containing additional inputs to likfn if needed, otherwise set to NULL

Details

These routines create lik and proc objects and call inneropt.

Value

A list with elements

coefs

Optimized coefficients at pars

lik

The lik object generated

proc

The proc item generated

res

The result of the optimization method

data

The data used in doing the fitting.

times

The vector of times at which the observations were made

See Also

inneropt, LS.setup, multinorm.setup, SplineCoefsErr

Examples

###############################
####   Data             #######
###############################

data(FhNdata)

###############################
####  Basis Object      #######
###############################

knots = seq(0,20,0.2)
norder = 3
nbasis = length(knots) + norder - 2
range = c(0,20)

bbasis = create.bspline.basis(range=range(FhNtimes),nbasis=nbasis,
	norder=norder,breaks=knots)


#### Start from pre-estimated values to speed up optimization

data(FhNest)

spars = FhNestPars
coefs = FhNestCoefs

lambda = 10000

res1 = Smooth.LS(make.fhn(),data=FhNdata,times=FhNtimes,pars=spars,coefs=coefs,
  basisvals=bbasis,lambda=lambda,in.meth='nlminb')


## Not run: 
# Henon system

hpars = c(1.4,0.3)              # Parameters
t = 1:200

x = c(-1,1)                     # Create some dataa
X = matrix(0,200+20,2)
X[1,] = x

for(i in 2:(200+20)){ X[i,] = make.Henon()$ode(i,X[i-1,],hpars,NULL) }

X = X[20+1:200,]

Y = X + 0.05*matrix(rnorm(200*2),200,2)

basisvals = diag(rep(1,200))    # Basis is just identiy
coefs = matrix(0,200,2)


# For sum of squared errors

lambda = 10000

res1	= Smooth.LS(make.Henon(),data=Y,times=t,pars=hpars,coefs,basisvals=basisvals,
  lambda=lambda,in.meth='nlminb',discrete=TRUE)

## End(Not run)


## Not run: 
# For multinormal transitions

var = c(1,0.01)

res2 = Smooth.multinorm(make.Henon(),data=Y,t,pars=hpars,coefs,basisvals=NULL,
  var=var,in.meth='nlminb',discrete=TRUE)

## End(Not run)