Permalink
Browse files

Compact implementation of sdv_multi example

  • Loading branch information...
1 parent 36802f2 commit 4e7c2adeab1d60c7113370060f562cadd970c065 @kaskr committed Feb 10, 2014
@@ -0,0 +1,34 @@
+require(TMB)
+
+# Read data
+source("sdv_multi_data.R")
+
+library(TMB)
+compile("sdv_multi_compact.cpp")
+dyn.load("sdv_multi_compact.so")
+obj <- MakeADFun(data=
+ list(n=n,p=p,y=t(y)),
+ parameters=list(
+ phi=rep(0.97,p),
+ log_sigma=rep(-1.7,p),
+ mu_x=rep(-0.5,p),
+ off_diag_x=rep(0.0,p),
+ h=t(matrix(0.0,nrow=n,ncol=p))
+ ),
+ random=c("^h")
+# map=list(off_diag_x=factor(rep(NA,3)))
+ )
+#obj$control <- list(trace=1,parscale=rep(1,12)*1e-2,REPORT=1,reltol=1e-12,maxit=100)
+
+obj$hessian <- F
+#ttt=obj$fn()
+##Rprof();opt <- do.call("optim",obj);Rprof(NULL)
+newtonOption(smartsearch=TRUE)
+##print(system.time(opt <- do.call("optim",obj)))
+ttt=system.time(opt<-nlminb(obj$par,obj$fn,obj$gr,
+ lower=c(rep(-.99,3),rep(-3.0,3),rep(-3.0,3),rep(-5.0,3)),
+ upper=c(rep(.99,3),rep(3.0,3),rep(3.0,3),rep(5.0,3))))
+##obj$gr()
+#c(phi1,phi2)
+#f(opt$par)
+rep <- sdreport(obj)
@@ -0,0 +1,33 @@
+// Compact version of sdv_multi
+#include <TMB.hpp>
+
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+ DATA_INTEGER(n);
+ DATA_INTEGER(p);
+ DATA_ARRAY(y);
+ PARAMETER_VECTOR(phi);
+ PARAMETER_VECTOR(log_sigma);
+ PARAMETER_VECTOR(mu_x);
+ PARAMETER_VECTOR(off_diag_x);
+ PARAMETER_ARRAY(h);
+
+ using namespace density;
+
+ Type g=0;
+ vector<Type> sigma=exp(log_sigma);
+ array<Type> ht=h.transpose(); // For row access
+ vector<Type> sigma_init=sigma/sqrt(Type(1.0)-phi*phi); // Initial sd of AR1
+ for(int j=0;j<p;j++)g += SCALE(AR1(phi(j)),sigma_init(j))(ht.col(j));
+
+ // Likelihood contribution: observations
+ UNSTRUCTURED_CORR_t<Type> neg_log_density(off_diag_x);
+ for(int i=0;i<n;i++)
+ {
+ vector<Type> sigma_y = exp(Type(0.5)*(mu_x + vector<Type>(h.col(i))));
+ g += VECSCALE(neg_log_density,sigma_y)(vector<Type>(y.col(i)));
+ }
+
+ return g;
+}
Binary file not shown.

0 comments on commit 4e7c2ad

Please sign in to comment.