Commit 40d92919 authored by kimwhoriskey's avatar kimwhoriskey

better way of adding in the gps errors

parent 7dbf4836
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR obj
using namespace density;
template<class Type>
Type dcrwSSMargos(objective_function<Type> * obj) {
// data
// names of each track
vector<std::string> tracknames = StringList("tracknames",obj);
int ntracks = tracknames.size();
vector<track<Type> > tracks(ntracks); // vector< kind of element> name(size)
// tracks
for(int i = 0; i < ntracks; i++){
track<Type> indtrack(tracknames[i],obj);
tracks[i] = indtrack;
}
REPORT(ntracks);
// parameters
// Input parameters - i.e. parameters to estimate.
PARAMETER_VECTOR(working_theta);
vector<Type> theta(2);
theta(0) = 2*M_PI/(1.0+exp(-working_theta(0))) - M_PI; // -pi to pi
theta(1) = 2*M_PI/(1.0+exp(-working_theta(1))); // 0 to 2pi
ADREPORT(theta);
PARAMETER_VECTOR(working_gamma);
vector<Type> gamma(2);
gamma(0) = 1.0/(1.0+exp(-working_gamma(0)));
gamma(1) = gamma(0)/(1.0+exp(-working_gamma(1)));
ADREPORT(gamma);
PARAMETER(working_sigma_lon);
Type sigma_lon = exp(working_sigma_lon);
ADREPORT(sigma_lon);
PARAMETER(working_sigma_lat);
Type sigma_lat = exp(working_sigma_lat);
ADREPORT(sigma_lat);
// Variance-covariance matrix for the process equation.
matrix<Type> Sigma(2,2);
Sigma << sigma_lon*sigma_lon, 0.0,
0.0, sigma_lat*sigma_lat;
PARAMETER_MATRIX(working_A); //matrix of switching probabilities
int m = working_A.rows(); // number of states
matrix<Type> A(m, m);
for(int i = 0; i < m; ++i){
for(int j = 0; j < (m-1); ++j){
A(i,j) = exp(working_A(i,j)); //ensures the probabilities we estimate are >0
}
A(i,m-1) = 1.0;
A.row(i) *= 1.0/A.row(i).sum();
}
ADREPORT(A);
// system of equations stat dist
matrix<Type> system(m, m);
for(int i=0; i<m; i++) {
for(int j=0; j<m; j++) {
if( i ==j ) {
system(i,j) = 1 - A(j,i) + 1;
} else {
system(i,j) = 0 - A(j,i) + 1;
}
}
}
// vector of ones for stat dist
vector<Type> ones(m);
for(int i=0; i<m; i++) {
ones(i) = 1.0;
}
// take the inverse of the system matrix and then multiply by the vector of ones
matrix<Type> sys_inverse = system.inverse();
vector<Type> delta = sys_inverse*ones;
PARAMETER(working_psi); //working value for measurement error
Type psi = exp(working_psi);
ADREPORT(psi);
// likelihood calculations
//Calculate logLikelihood
Type nll = 0.0;
// behaviour equation
vector<Type> behav(ntracks);
Type tmpb = 0.0;
for(int i=0; i < ntracks; ++i){
tmpb = behaviour(tracks[i].b, delta, A);
behav(i) = tmpb;
}
REPORT(behav);
nll -= behav.sum();
// movement equation
Type tmpp = 0.0;
vector<Type> move(ntracks);
for(int i=0; i < ntracks; ++i){
tmpp = movement(tracks[i].x, tracks[i].b, gamma, theta, Sigma);
move(i) = tmpp;
}
REPORT(move);
nll += move.sum();
// measurement equation
Type tmpm = 0.0;
vector<Type> meas(ntracks);
for(int i=0; i < ntracks; ++i){
tmpm = measurement_argos(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi);
meas(i) = tmpm;
}
REPORT(meas);
nll += meas.sum();
for(int i=0; i<ntracks; ++i) tracks[i].Rep();
REPORT(nll);
return nll;
}
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR this
......@@ -5,7 +5,7 @@ using namespace density;
template<class Type>
Type dcrwSSM(objective_function<Type> * obj) {
Type dcrwSSMcombo(objective_function<Type> * obj) {
// data
......@@ -82,15 +82,15 @@ Type dcrwSSM(objective_function<Type> * obj) {
matrix<Type> sys_inverse = system.inverse();
vector<Type> delta = sys_inverse*ones;
// PARAMETER(working_psi); //working value for measurement error
// Type psi = exp(working_psi);
// ADREPORT(psi);
//
// PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error
// vector<Type> gps_err = exp(working_gps_err);
// matrix<Type> gpsSigma(2,2);
// gpsSigma << gps_err(0)*gps_err(0), 0.0,
// 0.0, gps_err(1)*gps_err(1);
PARAMETER(working_psi); //working value for measurement error
Type psi = exp(working_psi);
ADREPORT(psi);
PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error
vector<Type> gps_err = exp(working_gps_err);
matrix<Type> gpsSigma(2,2);
gpsSigma << gps_err(0)*gps_err(0), 0.0,
0.0, gps_err(1)*gps_err(1);
......@@ -124,9 +124,13 @@ Type dcrwSSM(objective_function<Type> * obj) {
// measurement equation
Type tmpm = 0.0;
vector<Type> meas(ntracks);
// for(int i=0; i < ntracks; ++i){
// int measdatatype=tracks[i].datatype(0);
// tmpm = measurement(obj, tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, measdatatype);
// meas(i) = tmpm;
// }
for(int i=0; i < ntracks; ++i){
int measdatatype=tracks[i].datatype(0);
tmpm = measurement(obj, tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, measdatatype);
tmpm = measurement_combo(tracks[i].y, tracks[i].x, tracks[i].kind, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi, gpsSigma);
meas(i) = tmpm;
}
// for(int i=0; i < ntracks; ++i){
......
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR obj
using namespace density;
template<class Type>
Type dcrwSSMgps(objective_function<Type> * obj) {
// data
// names of each track
vector<std::string> tracknames = StringList("tracknames",obj);
int ntracks = tracknames.size();
vector<track<Type> > tracks(ntracks); // vector< kind of element> name(size)
// tracks
for(int i = 0; i < ntracks; i++){
track<Type> indtrack(tracknames[i],obj);
tracks[i] = indtrack;
}
REPORT(ntracks);
// parameters
// Input parameters - i.e. parameters to estimate.
PARAMETER_VECTOR(working_theta);
vector<Type> theta(2);
theta(0) = 2*M_PI/(1.0+exp(-working_theta(0))) - M_PI; // -pi to pi
theta(1) = 2*M_PI/(1.0+exp(-working_theta(1))); // 0 to 2pi
ADREPORT(theta);
PARAMETER_VECTOR(working_gamma);
vector<Type> gamma(2);
gamma(0) = 1.0/(1.0+exp(-working_gamma(0)));
gamma(1) = gamma(0)/(1.0+exp(-working_gamma(1)));
ADREPORT(gamma);
PARAMETER(working_sigma_lon);
Type sigma_lon = exp(working_sigma_lon);
ADREPORT(sigma_lon);
PARAMETER(working_sigma_lat);
Type sigma_lat = exp(working_sigma_lat);
ADREPORT(sigma_lat);
// Variance-covariance matrix for the process equation.
matrix<Type> Sigma(2,2);
Sigma << sigma_lon*sigma_lon, 0.0,
0.0, sigma_lat*sigma_lat;
PARAMETER_MATRIX(working_A); //matrix of switching probabilities
int m = working_A.rows(); // number of states
matrix<Type> A(m, m);
for(int i = 0; i < m; ++i){
for(int j = 0; j < (m-1); ++j){
A(i,j) = exp(working_A(i,j)); //ensures the probabilities we estimate are >0
}
A(i,m-1) = 1.0;
A.row(i) *= 1.0/A.row(i).sum();
}
ADREPORT(A);
// system of equations stat dist
matrix<Type> system(m, m);
for(int i=0; i<m; i++) {
for(int j=0; j<m; j++) {
if( i ==j ) {
system(i,j) = 1 - A(j,i) + 1;
} else {
system(i,j) = 0 - A(j,i) + 1;
}
}
}
// vector of ones for stat dist
vector<Type> ones(m);
for(int i=0; i<m; i++) {
ones(i) = 1.0;
}
// take the inverse of the system matrix and then multiply by the vector of ones
matrix<Type> sys_inverse = system.inverse();
vector<Type> delta = sys_inverse*ones;
PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error
vector<Type> gps_err = exp(working_gps_err);
matrix<Type> gpsSigma(2,2);
gpsSigma << gps_err(0)*gps_err(0), 0.0,
0.0, gps_err(1)*gps_err(1);
// likelihood calculations
//Calculate logLikelihood
Type nll = 0.0;
// behaviour equation
vector<Type> behav(ntracks);
Type tmpb = 0.0;
for(int i=0; i < ntracks; ++i){
tmpb = behaviour(tracks[i].b, delta, A);
behav(i) = tmpb;
}
REPORT(behav);
nll -= behav.sum();
// movement equation
Type tmpp = 0.0;
vector<Type> move(ntracks);
for(int i=0; i < ntracks; ++i){
tmpp = movement(tracks[i].x, tracks[i].b, gamma, theta, Sigma);
move(i) = tmpp;
}
REPORT(move);
nll += move.sum();
// measurement equation
Type tmpm = 0.0;
vector<Type> meas(ntracks);
for(int i=0; i < ntracks; ++i){
tmpm = measurement_gps(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, gpsSigma);
meas(i) = tmpm;
}
REPORT(meas);
nll += meas.sum();
for(int i=0; i<ntracks; ++i) tracks[i].Rep();
REPORT(nll);
return nll;
}
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR this
......@@ -43,29 +43,78 @@ Type movement(matrix<Type> x, vector<int> b, vector<Type> gamma, vector<Type> th
return rslt;
}
// // measurement process
// template<class Type>
// Type measurement(array<Type> y, matrix<Type> x, vector<int> idx, vector<Type> jidx, matrix<Type> ae, Type psi){
//
// vector<Type> measnll(y.cols());
// vector<Type> xhat(2);
// vector<Type> tmp(2);
// vector<Type> tmp2(2);
//
// for(int i = 0; i < y.cols(); ++i){
// xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i));
// tmp = y.col(i);
// tmp2 = tmp-xhat;
// // nll += nll_meas(tmp-xhat); //MVNORM_t returns the negative log likelihood
// // taken from Marie's DCRW_Argos.cpp file
// // dt returns the log likelihood
// measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude
// measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude
// }
//
// Type rslt = measnll.sum();
// return rslt; // negative log-lik
// }
// measurement process
template<class Type>
Type measurement_argos(array<Type> y, matrix<Type> x, vector<int> idx, vector<Type> jidx, matrix<Type> ae, Type psi){
vector<Type> measnll(y.cols());
vector<Type> xhat(2);
vector<Type> tmp(2);
vector<Type> tmp2(2);
for(int i = 0; i < y.cols(); ++i){
xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i));
tmp = y.col(i);
tmp2 = tmp-xhat;
// nll += nll_meas(tmp-xhat); //MVNORM_t returns the negative log likelihood
// taken from Marie's DCRW_Argos.cpp file
// dt returns the log likelihood
measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude
measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude
}
Type rslt = measnll.sum();
return rslt; // negative log-lik
}
template<class Type>
Type measurement_gps(array<Type> y, matrix<Type> x, vector<int> idx, vector<Type> jidx, matrix<Type> ae, matrix<Type> gpsSigma){
vector<Type> measnll(y.cols());
vector<Type> xhat(2);
vector<Type> tmp(2);
vector<Type> tmp2(2);
density::MVNORM_t<Type> gps_nll(gpsSigma);
for(int i = 0; i < y.cols(); ++i){
xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i));
tmp = y.col(i);
tmp2 = tmp-xhat;
measnll(i) = gps_nll(tmp2); //MVNORM_t returns the negative log likelihood
}
Type rslt = measnll.sum();
return rslt; // negative log-lik
}
template<class Type>
Type measurement_combo(array<Type> y, matrix<Type> x, vector<int> kind, vector<int> idx, vector<Type> jidx, matrix<Type> ae, Type psi, matrix<Type> gpsSigma){
vector<Type> measnll(y.cols());
vector<Type> xhat(2);
vector<Type> tmp(2);
vector<Type> tmp2(2);
density::MVNORM_t<Type> gps_nll(gpsSigma);
for(int i = 0; i < y.cols(); ++i){
xhat = (1.0-jidx(i)) * x.col(idx(i)-1) + jidx(i)*x.col(idx(i));
tmp = y.col(i);
tmp2 = tmp-xhat;
if(kind(i) == 0){
measnll(i) = gps_nll(tmp2); //MVNORM_t returns the negative log likelihood
} else if(kind(i) == 1) {
// taken from Marie's DCRW_Argos.cpp file
// dt returns the log likelihood
measnll(i) = -(0.5*log(psi) - log(ae(i,0)) + dt(sqrt(psi)*tmp2(0)/ae(i,0),ae(i,1),true) ); // Longitude
measnll(i) -= (0.5*log(psi) - log(ae(i,2)) + dt(sqrt(psi)*tmp2(1)/ae(i,2),ae(i,3),true) ); // Latitude
}
}
Type rslt = measnll.sum();
return rslt; // negative log-lik
}
//
......
......@@ -5,13 +5,17 @@
#include "likelihoods.hpp"
#include "measlikelihood.hpp"
#include "dcrws_hmm.hpp"
#include "dcrws_ssm.hpp"
#include "dcrws_ssm_combo.hpp"
#include "dcrws_ssm_argos.hpp"
#include "dcrws_ssm_gps.hpp"
// map the data_integer to a character
enum model_choices {
dcrwhmm = 0,
dcrwssm = 1,
dcrwssmargos = 1,
dcrwssmgps = 2,
dcrwssmcombo = 3,
};
template<class Type>
......@@ -24,8 +28,14 @@ switch(model){
case dcrwhmm:
return dcrwHMM(this);
break;
case dcrwssm:
return dcrwSSM(this);
case dcrwssmargos:
return dcrwSSMargos(this);
break;
case dcrwssmgps:
return dcrwSSMgps(this);
break;
case dcrwssmcombo:
return dcrwSSMcombo(this);
break;
}
}
......@@ -9,7 +9,7 @@ struct track {
array<Type> y;
vector<int> b;
vector<int> kind;
vector<int> datatype;
// vector<int> datatype;
vector<int> idx;
vector<Type> jidx;
matrix<Type> ae;
......@@ -30,7 +30,7 @@ track<Type>::track(std::string name, objective_function<Type> *obj) { // objecti
y = tmbutils::asArray<Type>(getListElement(tracking_data, "y"));
b = asVector<int>(getListElement(tracking_data, "b"));
kind = asVector<int>(getListElement(tracking_data, "kind"));
datatype = asVector<int>(getListElement(tracking_data, "datatype")) ;
// datatype = asVector<int>(getListElement(tracking_data, "datatype")) ;
idx = asVector<int>(getListElement(tracking_data, "idx"));
jidx = asVector<Type>(getListElement(tracking_data, "jidx"));
ae = asMatrix<Type>(getListElement(tracking_data, "ae"));
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment