dcrws_ssm_argos.hpp 3.25 KB
Newer Older
kimwhoriskey's avatar
kimwhoriskey committed
1 2 3 4 5 6 7
#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR obj

using namespace density;


template<class Type>
8
Type dcrwSSMargos(objective_function<Type> * obj) {
kimwhoriskey's avatar
kimwhoriskey committed
9

kimwhoriskey's avatar
kimwhoriskey committed
10
  // data
kimwhoriskey's avatar
kimwhoriskey committed
11

kimwhoriskey's avatar
kimwhoriskey committed
12 13 14 15
  // 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)
kimwhoriskey's avatar
kimwhoriskey committed
16

kimwhoriskey's avatar
kimwhoriskey committed
17 18 19 20 21
  // tracks
  for(int i = 0; i < ntracks; i++){
    track<Type> indtrack(tracknames[i],obj);
    tracks[i] = indtrack;
  }
kimwhoriskey's avatar
kimwhoriskey committed
22 23 24



kimwhoriskey's avatar
kimwhoriskey committed
25 26 27 28
  REPORT(ntracks);


  // parameters
kimwhoriskey's avatar
kimwhoriskey committed
29

kimwhoriskey's avatar
kimwhoriskey committed
30 31
  // Input parameters - i.e. parameters to estimate.
  PARAMETER_VECTOR(working_theta);
kimwhoriskey's avatar
kimwhoriskey committed
32 33 34
  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
kimwhoriskey's avatar
kimwhoriskey committed
35
  ADREPORT(theta);
kimwhoriskey's avatar
kimwhoriskey committed
36

kimwhoriskey's avatar
kimwhoriskey committed
37
  PARAMETER_VECTOR(working_gamma);
kimwhoriskey's avatar
kimwhoriskey committed
38 39 40
  vector<Type> gamma(2);
  gamma(0) = 1.0/(1.0+exp(-working_gamma(0)));
  gamma(1) = gamma(0)/(1.0+exp(-working_gamma(1)));
kimwhoriskey's avatar
kimwhoriskey committed
41
  ADREPORT(gamma);
kimwhoriskey's avatar
kimwhoriskey committed
42

kimwhoriskey's avatar
kimwhoriskey committed
43
  PARAMETER(working_sigma_lon);
kimwhoriskey's avatar
kimwhoriskey committed
44
  Type sigma_lon = exp(working_sigma_lon);
kimwhoriskey's avatar
kimwhoriskey committed
45
  ADREPORT(sigma_lon);
kimwhoriskey's avatar
kimwhoriskey committed
46
  PARAMETER(working_sigma_lat);
kimwhoriskey's avatar
kimwhoriskey committed
47
  Type sigma_lat = exp(working_sigma_lat);
kimwhoriskey's avatar
kimwhoriskey committed
48
  ADREPORT(sigma_lat);
kimwhoriskey's avatar
kimwhoriskey committed
49 50 51 52
  // 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;
kimwhoriskey's avatar
kimwhoriskey committed
53

kimwhoriskey's avatar
kimwhoriskey committed
54
  PARAMETER_MATRIX(working_A); //matrix of switching probabilities
kimwhoriskey's avatar
kimwhoriskey committed
55 56 57 58 59 60 61 62 63
  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();
  }
kimwhoriskey's avatar
kimwhoriskey committed
64 65
  ADREPORT(A);
  // system of equations stat dist
kimwhoriskey's avatar
kimwhoriskey committed
66 67 68 69 70 71 72 73 74 75
  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;
          }
      }
  }
kimwhoriskey's avatar
kimwhoriskey committed
76
  // vector of ones for stat dist
kimwhoriskey's avatar
kimwhoriskey committed
77 78 79 80 81 82 83 84
  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;

85 86 87
  PARAMETER(working_psi); //working value for measurement error
  Type psi = exp(working_psi);
  ADREPORT(psi);
kimwhoriskey's avatar
kimwhoriskey committed
88 89 90 91 92





kimwhoriskey's avatar
kimwhoriskey committed
93
   // likelihood calculations
kimwhoriskey's avatar
kimwhoriskey committed
94

95 96
  //Calculate logLikelihood
  Type nll = 0.0;
kimwhoriskey's avatar
kimwhoriskey committed
97

98
  // behaviour equation
kimwhoriskey's avatar
kimwhoriskey committed
99 100 101 102 103 104
  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;
  }
105
  REPORT(behav);
kimwhoriskey's avatar
kimwhoriskey committed
106
  nll -= behav.sum();
kimwhoriskey's avatar
kimwhoriskey committed
107

108
  // movement equation
kimwhoriskey's avatar
kimwhoriskey committed
109 110 111 112 113 114 115 116
  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();
kimwhoriskey's avatar
kimwhoriskey committed
117

118
  // measurement equation
kimwhoriskey's avatar
kimwhoriskey committed
119 120 121
  Type tmpm = 0.0;
  vector<Type> meas(ntracks);
  for(int i=0; i < ntracks; ++i){
122
    tmpm = measurement_argos(tracks[i].y, tracks[i].x, tracks[i].idx, tracks[i].jidx, tracks[i].ae, psi);
kimwhoriskey's avatar
kimwhoriskey committed
123 124 125 126 127
    meas(i) = tmpm;
  }
  REPORT(meas);
  nll += meas.sum();

kimwhoriskey's avatar
kimwhoriskey committed
128
  for(int i=0; i<ntracks; ++i) tracks[i].Rep();
kimwhoriskey's avatar
kimwhoriskey committed
129

130 131
  REPORT(nll);
  return nll;
kimwhoriskey's avatar
kimwhoriskey committed
132 133 134 135 136

}

#undef TMB_OBJECTIVE_PTR
#define TMB_OBJECTIVE_PTR this