Commit 7dbf4836 authored by kimwhoriskey's avatar kimwhoriskey

unfinished incorporation of the gps data, switching tactics now

parent c54df361
......@@ -82,9 +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(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);
......@@ -119,9 +125,14 @@ Type dcrwSSM(objective_function<Type> * obj) {
Type tmpm = 0.0;
vector<Type> meas(ntracks);
for(int i=0; i < ntracks; ++i){
tmpm = measurement(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi);
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){
// tmpm = measurement(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi);
// meas(i) = tmpm;
// }
REPORT(meas);
nll += meas.sum();
......
......@@ -43,31 +43,123 @@ 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(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
// #undef TMB_OBJECTIVE_PTR
// #define TMB_OBJECTIVE_PTR obj
// template<class Type>
// Type measurement(array<Type> y, matrix<Type> x, vector<int> kind, vector<int> idx, vector<Type> jidx, matrix<Type> ae, int datatype){
//
// vector<Type> measnll(y.cols());
// vector<Type> xhat(2);
// vector<Type> tmp(2);
// vector<Type> tmp2(2);
// // density::MVNORM_t<Type> gps_nll(gpsSigma);
//
// if(datatype==0){
//
// PARAMETER(working_psi); //working value for measurement error
// Type psi = exp(working_psi);
// ADREPORT(psi);
//
// 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;
//
// // 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
//
// }
//
// } else if (datatype==1){
//
// PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error
// vector<Type> gps_err = exp(working_gps_err);
// ADREPORT(gps_err);
// matrix<Type> gpsSigma(2,2);
// gpsSigma << gps_err(0)*gps_err(0), 0.0,
// 0.0, gps_err(1)*gps_err(1);
// 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
//
// }
//
//
// } else if (datatype==2) {
//
// 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);
// ADREPORT(gps_err);
// matrix<Type> gpsSigma(2,2);
// gpsSigma << gps_err(0)*gps_err(0), 0.0,
// 0.0, gps_err(1)*gps_err(1);
// 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;
//
// // different measurement error dists
// // 0 for gps, 1 for argos LS, could do 2 for Argos ellipse
// 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
// }
// #undef TMB_OBJECTIVE_PTR
// #define TMB_OBJECTIVE_PTR this
......
// measurement process
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR obj
template<class Type>
Type measurement(objective_function<Type> * obj, array<Type> y, matrix<Type> x, vector<int> kind, vector<int> idx, vector<Type> jidx, matrix<Type> ae, int datatype){
vector<Type> measnll(y.cols());
vector<Type> xhat(2);
vector<Type> tmp(2);
vector<Type> tmp2(2);
// density::MVNORM_t<Type> gps_nll(gpsSigma);
if(datatype==0){
PARAMETER(working_psi); //working value for measurement error
Type psi = exp(working_psi);
ADREPORT(psi);
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;
// 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
}
} else if (datatype==1){
PARAMETER_VECTOR(working_gps_err); // working value for gps measurement error
vector<Type> gps_err = exp(working_gps_err);
ADREPORT(gps_err);
matrix<Type> gpsSigma(2,2);
gpsSigma << gps_err(0)*gps_err(0), 0.0,
0.0, gps_err(1)*gps_err(1);
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
}
} else if (datatype==2) {
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);
ADREPORT(gps_err);
matrix<Type> gpsSigma(2,2);
gpsSigma << gps_err(0)*gps_err(0), 0.0,
0.0, gps_err(1)*gps_err(1);
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;
// different measurement error dists
// 0 for gps, 1 for argos LS, could do 2 for Argos ellipse
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
}
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR this
......@@ -3,6 +3,7 @@
#include "extramacros.hpp"
#include "trackstruct.hpp"
#include "likelihoods.hpp"
#include "measlikelihood.hpp"
#include "dcrws_hmm.hpp"
#include "dcrws_ssm.hpp"
......
......@@ -8,6 +8,8 @@ struct track {
track(std::string name,objective_function<Type> *obj);
array<Type> y;
vector<int> b;
vector<int> kind;
vector<int> datatype;
vector<int> idx;
vector<Type> jidx;
matrix<Type> ae;
......@@ -27,6 +29,8 @@ track<Type>::track(std::string name, objective_function<Type> *obj) { // objecti
SEXP tracking_data = getListElement(obj -> data,nam);
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")) ;
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