-
Notifications
You must be signed in to change notification settings - Fork 80
/
mvrw_sparse.cpp
32 lines (30 loc) · 1.13 KB
/
mvrw_sparse.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
// Identical with random walk example. Utilizing sparse block structure so efficient when the number of states is high.
#include <TMB.hpp>
/* Parameter transform */
template <class Type>
Type f(Type x){return Type(2)/(Type(1) + exp(-Type(2) * x)) - Type(1);}
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_ARRAY(obs); /* timeSteps x stateDim */
PARAMETER_ARRAY(u); /* State */
PARAMETER(transf_rho);
PARAMETER_VECTOR(logsds);
PARAMETER_VECTOR(logsdObs);
int timeSteps=obs.dim[1];
int stateDim=obs.dim[0];
Type rho=f(transf_rho);
vector<Type> sds=exp(logsds);
vector<Type> sdObs=exp(logsdObs);
// Setup object for evaluating multivariate normal likelihood
using namespace density;
VECSCALE_t<AR1_t<N01<Type> > > neg_log_density=VECSCALE(AR1(rho),sds);
/* Define likelihood */
Type ans=0;
ans-=dnorm(vector<Type>(u.col(0)),Type(0),Type(1),1).sum();
for(int i=1;i<timeSteps;i++)
ans+=neg_log_density(u.col(i)-u.col(i-1)); // Process likelihood
for(int i=1;i<timeSteps;i++)
ans-=dnorm(vector<Type>(obs.col(i)),vector<Type>(u.col(i)),sdObs,1).sum(); // Data likelihood
return ans;
}