BMAnimalTrack: Bayesian melding bias correction and uncertainty...

Description Usage Arguments Details Value Note Author(s) References Examples

Description

This is a wrapper function of function nllh.BB.Phi_XY, zSearch, and postMar.BB.Eta. It first estimate the parameters σ^2_H, σ^2_D by numerical maximizing nllh.BB.Phi_XY based on the GPS observations and DR path at the GPS time points. Then the grid of numerical integration of searched by zSearch. Finally, the posterior mean and variance of η is computed by postMar.BB.Eta. The posterior mean can be used as the corrected path. The posterior credible intervals are also calculated.

Usage

1
BMAnimalTrack(dataList, controlList=BMControl())

Arguments

dataList

A list produced by as.dataList

controlList

A list of parameters controls the whole calculation. Details are in BMControl

Details

This function carries out the Bayesian inference for the following model proposed in Liu et al.(2014):

η(t), the animal's path (one dimensional, northing or easting), is assumed to a Brownian Bridge process with variance parameter σ^2_H.(s2H)

Y, GPS observations observated at a subset of the time points of the path. Conditioning on η(t), the GPS observations Y(t) are unbiased iid normal observations of

Y(t) = η(t) + ε(t)

, where the measurment errors ε(t) iid N(0, σ^2_G). The σ^2_G (s2G) can be obtained from the manual of the GPS tags and is an input to this model.

X, Dead-Reckoned path, is assumed to be

X(t) = η(t) + h(t) + ξ(t),

where h(t) = ∑_{q=1}^Q d_q(t)β_q is a parametric function. Users can specify the design matrix [d_q(t)] or take the default polynomial d_q(t) = t^{q-1} and specify the order of it. ξ(t) is a zero-mean Brownian motion process, whose variance parameter is σ^2_D(s2D).

Value

When controlList$returnParam is false. Returns a T * 4 matrix. T is the number of times points in the path and the four columns of this matrix are:

Mean

Posterior mean of η

Variance

Posterior variance of η

CI.lower

The lower bound of the Bayesian credible intervals

CI.upper

The upper bound of the Bayesian credible intervals

When controlList$returnParam is true. The empirical Bayesian estimates of σ^2_H, σ^2_D(s2H, s2D) will also be returned.

Note

The current version cannot consider multiple initial values for the numerical maximization nor pass control parameters to the nlm used to perform the maximization.

Author(s)

Yang (Seagle) Liu <yang.liu@stat.ubc.ca>

References

Liu, Y., Battaile, B. C., Zidek, J. V., and Trites, A. (2014). Bayesian melding of the Dead-Reckoned path and gps measurements for an accurate and high-resolution path of marine mammals. arXiv preprint arXiv: 1411.6683.

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
set.seed(1)
#Generating data from our model
dlist <- dataSim(T=100, K=10, s2H=1, s2D=0.1, betaVec=c(1))
gpsObs <- dlist$Y
gpsTime <- dlist$Ytime
drPath <- dlist$X
#Produce the data list required by BMAnimalTrack
wlist <- as.dataList(drPath, gpsObs, gpsTime, 
            timeUnit=1, s2G=0.01, dUnit=1, betaOrder=1)
#Calculate the posterior of eta with our Bayesian Melding approach
etaMar <- BMAnimalTrack(wlist, BMControl(print=TRUE))

## Not run: 
## A real data example from package TrackReconstruction.
## Example from TrackReconstruction package.
## library(TrackReconstruction)
betas<-Standardize(1,1,-1,1,1,1,-57.8,68.76,-61.8,64.2,
    -70.16,58.08, -10.1,9.55,-9.75,9.72, -9.91,9.43)
#get declination and inclination data for study area
decinc<-c(10.228,65.918)
#data set with 6 associated GPS fixes in the "gpsdata" data set
data(rawdata)
#Perform the Dead-Reckoning of the raw accelerometer and 
# magnetometer data
DRoutput<-DeadReckoning(rawdata,betas,decinc,Hz=16,
            RmL=2,DepthHz=1,SpdCalc=3,MaxSpd=3.5)
#prepare GPS data
data(gpsdata02)
#matching time of the GPS and DR
gpsdata <- gpsdata02[gpsdata02$DateTime %in% DRoutput$DateTime, ]
gpsformat<-GPStable(gpsdata)
K <- nrow(gpsformat)
T <- nrow(gpsformat)
#Cut out the periods of DR path with the GPS
DRstart <- min(which(DRoutput$DateTime==gpsformat$DateTime[1]))
DRend <- max(which(DRoutput$DateTime==gpsformat$DateTime[K]))
#Thin the data (Original 16Hz, for now only working with 1Hz)
DRworking <- DRoutput[c(DRstart:DRend)[c(DRstart:DRend)%%16==1], ]

#Calculate the northing in km##
GPSnorthing=c(0, cumsum(gpsformat$DistanceKm[-1]*cos(gpsformat$BearingRad[-T])))
DRnorthing <- (DRworking$Ydim - DRworking$Ydim[1])/1000 
#Original unit of DR is in meters

#Data preparation for BM bias correction
ndl <- as.dataList(DRnorthing, GPSnorthing, 
    Ytime=gpsformat$DateTime, 
    Xtime=DRworking$DateTime, 
    s2G=0.0625, timeUnit=60, betaOrder=1)
#Bayesian Melding calculation.
nEtaMar <- BMAnimalTrack(ndl, BMControl(print=TRUE, returnParam=TRUE))
		
#Plots.
plot(ndl$XMx[,2], ndl$XMx[,1], type="l", col="blue", ylim=c(0, 2.5)) 
#uncorrected DR path
points(ndl$glist$Gtime, ndl$glist$Y, col="red", pch=16) #GPS points
lines(ndl$XMx[,2], nEtaMar$etaMar[,1], type="l") #Corrected path
lines(ndl$XMx[,2], nEtaMar$etaMar[,3], type="l", col="grey70") #Lower bound of CI
lines(ndl$XMx[,2], nEtaMar$etaMar[,4], type="l", col="grey70") #Upper bound of CI.

## End(Not run)

BayesianAnimalTracker documentation built on May 2, 2019, 5:39 a.m.