Permalink
Browse files

New example: theta logistic

  • Loading branch information...
Anders Nielsen
Anders Nielsen committed Feb 7, 2014
1 parent dc13562 commit 7734d66d8aecb6ce259fda9b8920290cbbeb3c93
@@ -0,0 +1,4 @@
+## Number of data points
+200
+## Data points
+ 3.023146 3.219844 3.654241 3.604655 3.542303 4.431661 4.291311 4.367869 4.37504 3.945371 4.291201 4.59584 4.302166 4.151933 4.428207 4.305366 4.548615 4.632054 4.784911 4.799618 4.926257 5.136269 5.148587 4.648586 5.29066 4.678948 5.080761 5.223667 5.261118 4.760974 4.789399 5.221183 5.207064 5.20803 5.182091 5.412485 5.815396 5.20875 5.787764 6.08957 5.900544 5.940576 6.214075 5.841425 6.420394 6.466935 6.737423 6.764397 6.806826 6.331149 6.779501 6.721639 6.741977 7.00933 6.658492 6.518948 6.295724 6.34795 6.683175 6.914437 6.275365 6.417207 6.826017 6.587109 6.505463 6.859157 6.678419 6.835149 6.888164 6.829884 6.623776 6.8927 6.609411 6.802324 6.888346 6.928283 7.121756 6.661365 6.7966 6.429909 6.591986 6.775834 6.491397 6.525561 6.956582 6.571743 6.62285 6.405347 6.524295 6.645468 6.148118 6.581335 6.248974 6.047809 6.51518 6.215392 6.994722 6.404407 6.997952 6.822027 6.735724 6.748609 6.308172 6.535844 6.370927 6.475951 6.292657 6.353783 6.538954 6.243774 6.493041 6.300049 6.568226 6.525136 6.751022 6.270681 6.833939 6.488327 6.235697 6.718966 6.614604 6.713407 6.620666 6.501483 6.384192 6.727767 6.769247 6.66819 6.311557 7.059169 6.690371 6.09881 6.803202 6.308663 6.563545 6.390233 6.59533 6.692878 6.682184 6.919849 6.900628 6.514584 6.882142 7.049829 7.076877 7.127657 6.595289 6.714436 6.956282 6.493041 6.764563 6.449598 6.522969 6.47359 6.720563 6.671342 6.198196 6.478438 6.918452 6.791604 6.964608 6.726904 6.490726 6.975198 7.013087 7.254823 6.792741 7.050987 6.951828 6.982148 7.106363 6.980717 7.327641 7.017561 6.512386 6.987826 6.643562 6.671664 7.217902 6.732841 6.509157 6.327128 6.901743 6.584902 6.719805 6.441681 6.742182 6.641601 6.875713 6.726997 7.130124 7.014675 6.922911 6.680622 7.072831 6.649716 7.01744 6.989642 6.710071 6.835827
@@ -0,0 +1,38 @@
+DATA_SECTION
+ init_int N
+ init_vector Y(1,N)
+
+PARAMETER_SECTION
+ init_number logr0
+ init_number logtheta
+ init_bounded_number logK(4.6,7.6)
+ init_number logQ
+ init_number logR
+
+ random_effects_vector X(1,N);
+
+ sdreport_number r0
+ sdreport_number theta
+ sdreport_number K
+ sdreport_number Q
+ sdreport_number R
+
+ objective_function_value jnll
+
+PROCEDURE_SECTION
+ for(int i=2; i<=N; ++i){
+ step(X(i-1),X(i),logr0,logK,logtheta,logQ); }
+
+ for(int i=1; i<=N; ++i){
+ obs(X(i),logR,i);}
+
+ r0=exp(logr0); theta=exp(logtheta); K=exp(logK); Q=exp(logQ); R=exp(logR);
+
+SEPARABLE_FUNCTION void step(const dvariable& x1, const dvariable& x2, const dvariable& logr0, const dvariable& logK, const dvariable& logtheta, const dvariable& logQ)
+ dvariable var=exp(logQ);
+ dvariable m=x1 + exp(logr0) * (1.0 - pow(exp(x1)/exp(logK),exp(logtheta)));
+ jnll+=0.5*(log(2.0*M_PI*var)+square(x2-m)/var);
+
+SEPARABLE_FUNCTION void obs(const dvariable& x, const dvariable& logR, int i)
+ dvariable var=exp(logR);
+ jnll+=0.5*(log(2.0*M_PI*var)+square(x-Y(i))/var);
View
@@ -0,0 +1,26 @@
+Y<-scan('thetalog.dat', skip=3, quiet=TRUE)
+
+library(TMB)
+compile("thetalog.cpp")
+dyn.load("thetalog.so")
+data <- list(Y=Y)
+parameters <- list(
+ X=data$Y*0,
+ logr0=0,
+ logtheta=0,
+ logK=6,
+ logQ=0,
+ logR=0
+ )
+newtonOption(smartsearch=FALSE)
+obj <- MakeADFun(data,parameters,random="^X",DLL="thetalog")
+obj$hessian <- TRUE
+#obj$control$reltol<-1e-12
+
+obj$fn()
+obj$gr()
+system.time(opt <- nlminb(obj$par,obj$fn,obj$gr))
+
+rep <- sdreport(obj)
+rep
+
View
@@ -0,0 +1,32 @@
+// Random walk with multivariate correlated increments and measurement noise.
+#include <TMB.hpp>
+
+/* Parameter transform */
+template <class Type>
+Type square(Type x){return x*x;}
+
+template<class Type>
+Type objective_function<Type>::operator() ()
+{
+ DATA_VECTOR(Y); /* timeSteps x stateDim */
+ PARAMETER_VECTOR(X); /* State */
+ PARAMETER(logr0);
+ PARAMETER(logtheta);
+ PARAMETER(logK);
+ PARAMETER(logQ)
+ PARAMETER(logR)
+
+ int timeSteps=Y.size();
+ using namespace density;
+ Type ans=0;
+ for(int i=1;i<timeSteps;i++){
+ Type var=exp(logQ);
+ Type m=X[i-1] + exp(logr0) * (1.0 - pow(exp(X[i-1])/exp(logK),exp(logtheta)));
+ ans+=0.5*(log(2.0*M_PI*var)+square(X[i]-m)/var);
+ }
+ for(int i=0;i<timeSteps;i++){
+ Type var=exp(logR);
+ ans+=0.5*(log(2.0*M_PI*var)+square(X[i]-Y[i])/var);
+ }
+ return ans;
+}
@@ -0,0 +1,4 @@
+## Number of data points
+200
+## Data points
+ 3.023146 3.219844 3.654241 3.604655 3.542303 4.431661 4.291311 4.367869 4.37504 3.945371 4.291201 4.59584 4.302166 4.151933 4.428207 4.305366 4.548615 4.632054 4.784911 4.799618 4.926257 5.136269 5.148587 4.648586 5.29066 4.678948 5.080761 5.223667 5.261118 4.760974 4.789399 5.221183 5.207064 5.20803 5.182091 5.412485 5.815396 5.20875 5.787764 6.08957 5.900544 5.940576 6.214075 5.841425 6.420394 6.466935 6.737423 6.764397 6.806826 6.331149 6.779501 6.721639 6.741977 7.00933 6.658492 6.518948 6.295724 6.34795 6.683175 6.914437 6.275365 6.417207 6.826017 6.587109 6.505463 6.859157 6.678419 6.835149 6.888164 6.829884 6.623776 6.8927 6.609411 6.802324 6.888346 6.928283 7.121756 6.661365 6.7966 6.429909 6.591986 6.775834 6.491397 6.525561 6.956582 6.571743 6.62285 6.405347 6.524295 6.645468 6.148118 6.581335 6.248974 6.047809 6.51518 6.215392 6.994722 6.404407 6.997952 6.822027 6.735724 6.748609 6.308172 6.535844 6.370927 6.475951 6.292657 6.353783 6.538954 6.243774 6.493041 6.300049 6.568226 6.525136 6.751022 6.270681 6.833939 6.488327 6.235697 6.718966 6.614604 6.713407 6.620666 6.501483 6.384192 6.727767 6.769247 6.66819 6.311557 7.059169 6.690371 6.09881 6.803202 6.308663 6.563545 6.390233 6.59533 6.692878 6.682184 6.919849 6.900628 6.514584 6.882142 7.049829 7.076877 7.127657 6.595289 6.714436 6.956282 6.493041 6.764563 6.449598 6.522969 6.47359 6.720563 6.671342 6.198196 6.478438 6.918452 6.791604 6.964608 6.726904 6.490726 6.975198 7.013087 7.254823 6.792741 7.050987 6.951828 6.982148 7.106363 6.980717 7.327641 7.017561 6.512386 6.987826 6.643562 6.671664 7.217902 6.732841 6.509157 6.327128 6.901743 6.584902 6.719805 6.441681 6.742182 6.641601 6.875713 6.726997 7.130124 7.014675 6.922911 6.680622 7.072831 6.649716 7.01744 6.989642 6.710071 6.835827

0 comments on commit 7734d66

Please sign in to comment.