# Simulation and interpolation of trajectories. `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. ```r 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. ## Returns - **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 ```r # 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) ```

sim_trajectory_joint function