# inst/MatlabCollocInfer/Rcode/SSElik.R In CollocInfer: Collocation Inference for Dynamic Systems

```############ squared measurement error ##############

#  This file defines functions for evaluating the total SSE for data
#  and its derivative values:

#    fn      = SSE
#    dfdx    = dSSE.dx
#    dfdy    = dSSE.dy
#    dfdp    = dSSE.dp
#    d2fdx2  = d2SSE.dx2
#    d2fdxdy = d2SSE.dxdy
#    d2fdy2  = d2SSE.dy2
#    d2fdxdp = d2SSE.dxdp
#    d2fdydp = d2SSE.dydp

# as well as function checkweights

################################################################################

make.SSElik <- function()
{

##############################################################################

SSE <- function(data,times,devals,pars,more)
{
#  Define the values of the process at observation times
fdevals = more\$fn(times,devals,pars,more\$more)
difs    = data - fdevals
#  if process is unobserved, data[i] is NaN, and so is difs[i]
difs[is.na(difs)] = 0
#  define weights
weights = checkweights(more\$weights,more\$whichobs,difs)
#  compute total weighted error sum of squares
f = apply( weights*difs^2, 1, sum )
return(f)
}

##############################################################################

dSSE.dx <- function(data,times,devals,pars,more)
{
fdevals = more\$fn(times,devals,pars,more\$more)
difs    = data-fdevals
difs[is.na(difs)] = 0
weights = checkweights(more\$weights,more\$whichobs,difs)
difs = weights*difs
dfdx = more\$dfdx(times,devals,pars,more\$more)
g = c()
for(i in 1:dim(dfdx)[3]){
g = cbind(g,apply(difs*dfdx[,,i],1,sum))
}
return(-2*g)
}

##############################################################################

dSSE.dy <- function(data,times,devals,pars,more)
{
fdevals = more\$fn(times,devals,pars,more\$more)
difs    = data-fdevals
difs[is.na(difs)] = 0
weights = checkweights(more\$weights,more\$whichobs,difs)
difs    = weights*difs
return(2*difs)
}

##############################################################################

dSSE.dp <- function(data,times,devals,pars,more)
{
fdevals = more\$fn(times,devals,pars,more\$more)
difs    = data-fdevals
difs[is.na(difs)] = 0
weights = checkweights(more\$weights,more\$whichobs,difs)
difs    = weights*difs
dfdp = more\$dfdp(times,devals,pars,more\$more)
g = c()
for(i in 1:dim(dfdp)[3]){
g  = cbind(g,apply(difs*dfdp[,,i],1,sum))
}
return(-2*g)
}

##############################################################################

d2SSE.dx2 <- function(data,times,devals,pars,more)
{
fdevals = more\$fn(times,devals,pars,more\$more)
difs    = data-fdevals
dfdx    = more\$dfdx(times,devals,pars,more\$more)
d2fdx2  = more\$d2fdx2(times,devals,pars,more\$more)
difs[is.na(difs)] = 0
weights = checkweights(more\$weights,more\$whichobs,difs)
difs   = weights*difs
H = array(0,c(dim(devals),dim(devals)[2]))
for(i in 1:dim(d2fdx2)[3]){
for(j in 1:dim(d2fdx2)[4]){
H[,i,j] = apply(-difs*d2fdx2[,,i,j] +  weights*dfdx[,,j]*dfdx[,,i],1,sum)
}
}

return(2*H)
}

##############################################################################

d2SSE.dxdy <- function(data,times,devals,pars,more)
{
dfdx   = more\$dfdx(times,devals,pars,more\$more)
weights = checkweights(more\$weights,more\$whichobs,dfdx[,,1])
weights[is.na(data)] = 0
for(i in 1:dim(dfdx)[3]){
dfdx[,,i] = weights*dfdx[,,i]
}
return(-2*aperm(dfdx,c(1,3,2)))
}

##############################################################################

d2SSE.dy2 <- function(data,times,devals,pars,more)
{
r = array(0,c(dim(data),dim(data)[2]))
ind = cbind(rep(1:dim(data)[1],dim(data)[2]),
kronecker(cbind(1:dim(data)[2],1:dim(data)[2]),rep(1,dim(data)[1])))
weights = checkweights(more\$weights,more\$whichobs,data)
r[ind]  = weights
return(2*r)
}

##############################################################################

d2SSE.dxdp <- function(data,times,devals,pars,more)
{
fdevals = more\$fn(times,devals,pars,more\$more)
difs = data-fdevals
difs[is.na(difs)] = 0
weights = checkweights(more\$weights,more\$whichobs,difs)
difs = weights*difs
dfdx = more\$dfdx(times,devals,pars,more\$more)
dfdp = more\$dfdp(times,devals,pars,more\$more)
d2fdxdp = more\$d2fdxdp(times,devals,pars,more\$more)
H = array(0,c(dim(devals),length(pars)))
for(i in 1:dim(d2fdxdp)[3]){
for(j in 1:dim(d2fdxdp)[4]){
H[,i,j] = apply(-difs*d2fdxdp[,,i,j] + weights*dfdx[,,i]*dfdp[,,j],1,sum)
}
}
return(2*H)
}

##############################################################################

d2SSE.dydp <- function(data,times,devals,pars,more)
{
dfdp = more\$dfdp(times,devals,pars,more\$more)
weights = checkweights(more\$weights,more\$whichobs,dfdp[,,1])
weights[is.na(data)] = 0
for(i in 1:dim(dfdp)[3]){
dfdp[,,i] = weights*dfdp[,,i]
}
return(-2*dfdp)
}

##############################################################################

return(
list(
fn      = SSE,
dfdx    = dSSE.dx,
dfdy    = dSSE.dy,
dfdp    = dSSE.dp,
d2fdx2  = d2SSE.dx2,
d2fdxdy = d2SSE.dxdy,
d2fdy2  = d2SSE.dy2,
d2fdxdp = d2SSE.dxdp,
d2fdydp = d2SSE.dydp
)
)
}

################################################################################

checkweights <- function(weights,whichrows,diffs){
if(is.null(whichrows)){ whichrows = 1:nrow(diffs) }

if(is.null(weights)){ return(matrix(1,dim(diffs))) }
else if( length(weights) == ncol(diffs)){ return( matrix(weights,nrow(diffs),ncol(diffs),byrow=TRUE)) }
else if( length(weights[whichrows]) == nrow(diffs)){ return( matrix(weights,nrow(diffs),ncol(diffs),byrow=FALSE)) }
else if( prod( dim(weights[whichrows,]) == dim(diffs)) ){ return(weights[whichrows,]) }
else if( ncol(weights[whichrows,]) > ncol(diffs) & nrow(weights) > nrow(diffs) ){
warning('Dimension of weights does not match that of data')
return( weights[whichrows[1:nrow(diffs)],1:ncol(diffs)] )
}
else{
stop('Dimension of weights does not match that of data')
}

}
```

## Try the CollocInfer package in your browser

Any scripts or data that you put into this service are public.

CollocInfer documentation built on May 2, 2019, 4:03 a.m.