sim_trajectory_joint {animalEKF}

Simulation and interpolation of trajectories.

Description

sim_trajectory_joint simulates regular-step trajectories under correlated random walk (CRW). interp_trajectory_joint interpolates regular steps to irregular ones drawn from a log-normal distribution.

Usage

sim_trajectory_joint(area_map=NULL, centroids=matrix(c(0,0), ncol=2), 
                     transition_matrices=list(matrix(c(10,3,2,9), 
                     ncol=2, byrow=TRUE)), 
                     mu0_pars=list(alpha=c(-4 ,-1.6), beta=c(0,0)), 
                     var0_pars=list(alpha=c(1.6,0.16), beta=c(2,0.5)), 
                     N=100, nstates=2, reg_dt=120, gen_irreg=TRUE, 
                     one_d=FALSE, dt_lnorm_mu=log(120), dt_lnorm_sd=1,
                     dt_vals=NULL, starting_polygon=area_map,
                     nsharks=1, interact=FALSE, 
                     interact_pars=list(interacting_sharks=c(1:nsharks),
                     time_radius=60*30, spat_radius=200, min_num_neibs=10,
                     eta_mu=c(2,1), rho_sd=c(0.75, 0.75)),
                     time_dep_trans=FALSE, trans_alpha=c(1, 1.5)) 

interp_trajectory_joint(d, nstates, one_d, dt_lnorm_mu=5, dt_lnorm_sd=1, 
                        dt_vals=NULL, centroids=matrix(c(0,0), ncol=2))

Arguments

area_map

Shapefile within which the observations are located (optional). Should be the output of applying sf::st_geometry on an object of class sf. If input is NULL, a default rectangular shapefile is created.

centroids

Matrix with two columns specifying the centroids of regions. The number of rows specifies the number of regions.

transition_matrices

A list of 2x2 matrices specifying the Dirichlet parameters for behavior transition probabilities. The list is replicated so it's the length of the number of regions. If nstates=1 then these are not used since there is only one behavior.

mu0_pars

List of mean values of alpha (=log-speed if 2-D, and velocity if 1-D) and beta (turn angle, ignored for 1-D) for one or two behavioral states.

var0_pars

List of variances of alpha and beta distributions (see mu0_pars).

N

Number of regular steps to simulate.

nstates

Number of behavioral states. For now restricted to a maximum of 2.

reg_dt

Length in seconds of each regular interval.

gen_irreg

Logical. If TRUE, then use interp_trajectory_joint to make irregular steps.

one_d

Logical. If TRUE, then simulation occurs on 1-D line, if FALSE (the default) it is 2-D.

dt_lnorm_mu

Mean parameter mu of the log-normal distribution to draw time step lengths.

dt_lnorm_sd

Standard deviation parameter sigma of the log-normal distribution to draw time step lengths.

starting_polygon

Polygon to draw starting coordinates in. This helps if you want the trajectories to start around the same area.

nsharks

Number of sharks to simulate trajectories for. If nsharks>1, then joint effects may take place.

interact

Logical. If TRUE, simulate interaction parameters of neighborhood (either 1-D or 2-D). If nstates=1, automatically set to FALSE.

interact_pars

List of interaction priors: 1) interacting_sharks means which of the sharks 1...nsharks are to use interaction parameters; 2) time_radius is the time in seconds, and 3) spat_radius is the spatial radius is meters to consider for spatial neighbors; 4) min_num_neibs is the minimum number of time and spatial radius observations that need to exist to constitute a neighborhood; 5) eta_mu is the vector of mean value for the interaction parameter eta; rho_sd is the vector of standard deviations of the interaction multiplier rho.

time_dep_trans

Logical. If TRUE, state transition matrices are time-dependent meaning that probability depends on the number of steps a shark has remained in the current state.

trans_alpha

If time_dep_trans=TRUE, the transition alpha parameters for the Dirichlet distribution for drawing behaviors.

d

Input for interp_trajectory_joint. An array, usually output by sim_trajectory_joint, of regular-step trajectories.

dt_vals

An optional vector of time difference values. By default is NULL, meaning time gaps will be generated by dt_lnorm_mu and dt_lnorm_sd, but supplying a vector to dt_vals lets the user specify the time gaps rather than having them be randomly generated.

Value

d

Array of regular-step trajectory locations.

d_ds

Object d in format data.frame.

di

If gen_irreg==TRUE, is the non-constant step length locations.

Author(s)

Samuel Ackerman

Examples

	

# read shapefile and convert into a Polygon
bolsachica <- sf::st_read(system.file("shapes/FTB_lines.shp", package="animalEKF")[1])
bolsachica <- sf::st_polygonize(bolsachica)

island <- sf::st_read(system.file("shapes/FTB_island.shp", package="animalEKF")[1])
# the actual available room for movement is the area in the water, subtracting the island inside
bolsachica <- sf::st_difference(sf::st_geometry(bolsachica), sf::st_geometry(island))


# sample 5 points approximately equally spaced within the shapefile, as region centroids
regions <- sf::st_sample(x=sf::st_geometry(bolsachica), size=5, type="regular", exact=TRUE)
# extract the coordinates
regions <- as.data.frame(sf::st_coordinates(regions)[, c("X","Y")])

#define Voronoi tessellation tile in which to start shark paths
vortess <- deldir::deldir(x=regions[,"X"], y=regions[,"Y"], wlines="tess", 
						  plotit=FALSE, suppressMsge=TRUE)
# convert these to a set of Polygons, and choose one of them as the starting polygon
vtiles <- sf::st_as_sf(tess2spat(vortess))
sf::st_crs(vtiles) <- sf::st_crs(bolsachica)
# extract only the 3rd tile
# note, want to have the simulation paths spread out, so a given draw may result in
# cramped and thus hard to estimate paths

starting_polygon <- sf::st_sfc(vtiles[[1]][[3]], crs=sf::st_crs(vtiles))
starting_polygon <- sf::st_intersection(sf::st_geometry(bolsachica),
						sf::st_geometry(starting_polygon))

#define list of transition matrices between behaviors

tmat_list <- list(matrix(c(8, 2, 2, 4), ncol=2, byrow=TRUE),
                  matrix(c(1.5*5, 1.5*1, 3, 3), ncol=2, byrow=TRUE),
                  matrix(c(7, 1, 1, 7), ncol=2, byrow=TRUE))

#generate 4-shark simulated trajectory with 100 regular steps of length 120 seconds.  
#Sharks 3 and 4 will be interacting with the others, but 1 and 2 will not.
			  
nsharks <- 4	
		
#simulate trajectory
#setting gen_irreg=TRUE generates an irregular trajectory from the regular-step one
#with the log-normal specified in dt_lnorm_mu and dt_lnorm_sd
#sim_4sharks$di would contain the irregular dataset
#otherwise, say you wanted to try different interpolations, you can use the same regular
#step from sim_trajectory_joint and then interpolate separately with interp_trajectory_joint.

#make simulated trajectories all start in the same area so they will be close enough to be
#interacting, for the purposes of this exercise
#note that the simulation may time out trying to draw points in this starting polygon that end 
#up in the shapefile boundary

nsteps_sim <- 100
reg_dt <- 120



sim_4sharks <- sim_trajectory_joint(area_map=sf::st_geometry(bolsachica), centroids=regions, 
                                    transition_matrices=tmat_list, nsharks=nsharks,
                                    mu0_pars=list(alpha=c(-4 ,-1.6), beta=c(0,0)), 
                                    var0_pars=list(alpha=c(1,0.25), beta=c(1,.25)),
                                    N=nsteps_sim, nstates=2, reg_dt=reg_dt,
                                    gen_irreg=FALSE, one_d=FALSE,
                                    starting_polygon=starting_polygon, interact=TRUE, 
                                    interact_pars=list(interacting_sharks=c(3:4),
                                    time_radius=60*30, spat_radius=150,
                                    min_num_neibs=10,
                                    eta_mu=c(2,1), rho_sd=c(0.75, 0.75)),
                                    time_dep_trans=FALSE, 
                                    dt_lnorm_mu=log(120), dt_lnorm_sd=0.4) 




#plot trajectories

shark_names <- dimnames(sim_4sharks$d)[[ 3 ]]
shark_colors <- 2:5
names(shark_colors) <- shark_names


sp::plot(bolsachica, main="Full trajectories")
deldir::plot.deldir(vortess, wlines="tess", add=TRUE)

for (ss in shark_names) {
    lines(sim_4sharks$d[,c("X","Y"), ss], col=shark_colors[ss])
}
		   
#now interpolate to uneven steps with lognormal mean log(120) (so they are on 
#average the same as the regular steps and sd=0.4
#d is the regular step, di is irregular

#if want to interpolate separately.  Otherwise just set gen_irreg=TRUE above
#this is so you can interpolate a dataset not generated by sim_trajectory_joint
#if gen_irreg=TRUE in sim_trajectory_joint,
#interp_ds will be returned as the 'di' object
				   
interp_ds <- interp_trajectory_joint(d=sim_4sharks$d, nstates=2, 
                                     one_d=FALSE, 
                                     dt_lnorm_mu=log(reg_dt),
                                     dt_lnorm_sd=0.4,
                                     centroids=regions)	



#now plot observed ones, may differ
sp::plot(bolsachica, main="Observed trajectories")
deldir::plot.deldir(vortess, wlines="tess", add=TRUE)

for (ss in shark_names) {
    lines(interp_ds[ interp_ds$tag == ss ,c("X","Y")], col=shark_colors[ss])
}


	
		
#try to recover EKF with steps at the original 120 seconds
#use the original simulated transition and foraging probabilities for comparison		
		
		
#intial values for some parameters
tau_pars_init <- c(8, 14, 10,1) #2
sigma_pars_init <- c(5, 8, 8, 3)

#measurement error

bmat <- matrix(c(1, -0.3, -0.3, 1), ncol=2)
Errvar_init1 <-5*20*bmat
Errvar_init2 <- 15*20*bmat

#particle error
Particle_err_init <- 0.5*20*bmat
		
# only estimate movement on first 5 steps
# for better results, npart should be set higher, like 150 or more
nsteps_estimate <- 5
npart <- 15		
		
		
#again, if you use gen_irreg=TRUE in sim_trajectory_joint,
#the input 'd' argument should be sim_4sharks$di or interp_ds

#NOTE: user should set output_plot=TRUE to see PDF,
#for purposes of package testing we set it to FALSE 
# if FALSE, plots will still appear in the console


ekf_interp_mod <- EKF_interp_joint(d=interp_ds, npart=npart, 
                                   area_map=bolsachica,
                                   state_favor=c(1,2),
                                   centroids=regions, 
                                   sigma_pars=sigma_pars_init,
                                   tau_pars=tau_pars_init,
                                   Errvar0=list(Errvar_init1, Errvar_init2),
                                   Particle_errvar0=Particle_err_init,
                                   mu0_pars=list(alpha=c(-4 ,-1.3), beta=c(0,0)), 
                                   truncate=TRUE, 
                                   neff_sample=0.75, dirichlet_init=c(8,2,2,4), 
                                   smoothing=TRUE, fix_smoothed_behaviors=FALSE, 
                                   time_dep_trans=FALSE, resamp_full_hist=FALSE, 
                                   nstates=2, reg_dt=reg_dt, interact=TRUE, 
                                   maxStep=nsteps_estimate, update_eachstep=TRUE, 
                                   compare_with_known=TRUE, 
                                   known_trans_prob=sim_4sharks$true_transition_prob, 
                                   known_foraging_prob=sim_4sharks$true_foraging_prob, 
                                   known_regular_step_ds=sim_4sharks$d_ds,
                                   output_plot=FALSE)		
				   

#simulate one-dimensional movement for 1 robot (shark)
#here we use gen_irreg=TRUE instead of generating a separate interpolation object

one_d <- sim_trajectory_joint(centroids=NULL, N=nsteps_sim,
                              mu0_pars=list(alpha=c(4, 9)),
                              var0_pars=list(alpha=c(1, 1)),
                              transition_matrices=tmat_list[[ 1 ]], nstates=2,
                              reg_dt=reg_dt, gen_irreg=TRUE, one_d=TRUE,
                              dt_lnorm_mu=log(120), dt_lnorm_sd=0.55)




#measurement error
bmat <- matrix(1)
Errvar_init1 <-1*bmat
Errvar_init2 <-3*bmat


#particle error
Particle_err_init <- 0.1*bmat



ekf_1d <- EKF_1d_interp_joint(d=one_d$di, npart=npart,  maxStep=nsteps_estimate, 
                              state_favor=c(1,1), nstates=2, lowvarsample=TRUE, 
                              neff_sample=1, time_dep_trans=FALSE, reg_dt=reg_dt, 
                              max_int_wo_obs=15, resamp_full_hist=FALSE, 
                              alpha0_pars=list(mu0=c(4, 9), V0=c(0.25, 0.25)),
                              sigma_pars=sigma_pars_init,
                              Errvar0=list(Errvar_init1, Errvar_init2),
                              Particle_errvar0=Particle_err_init,
                              compare_with_known=TRUE, 
                              known_trans_prob=one_d$true_transition_prob, 
                              known_foraging_prob=one_d$true_foraging_prob, 
                              known_regular_step_ds=one_d$d_ds, update_eachstep=TRUE,
                              smoothing=TRUE, output_plot=FALSE)