Set-up (1) Reading data
data <- read.table("linreg.dat", header=TRUE)
data
## x y
## 1 -1 1.4
## 2 0 4.7
## 3 1 5.1
## 4 2 8.3
## 5 3 9.0
## 6 4 14.5
## 7 5 14.0
## 8 6 13.4
## 9 7 19.2
## 10 8 18.0
Set-up (2) Writing cpp file as “linreg.cpp”"
#include <TMB.hpp>
template<class Type>
Type objective_function<Type>::operator() ()
{
DATA_VECTOR(x);
DATA_VECTOR(y);
int n = y.size();
PARAMETER(b0);
PARAMETER(b1);
PARAMETER(logSigma);
vector<Type> yfit(n);
Type neglogL = 0.0;
yfit = b0 + b1*x;
neglogL = -sum(dnorm(y, yfit, exp(logSigma), true));
return neglogL;
}
Set-up (3) Compiling “cpp” file and call “Dynamic link”
parameters <- list(b0=0, b1=0, logSigma=0)
require(TMB)
## Loading required package: TMB
compile("linreg.cpp")
## [1] 0
dyn.load(dynlib("linreg"))
Minimizing an objective function defined
- “MakeADFun”: Constructing objective functions with derivatives based on a compiled C++ template
- “nlminb”: Minimization
model <- MakeADFun(data, parameters, DLL="linreg")
fit <- nlminb(model$par, model$fn, model$gr)
## outer mgc: 1468.4
## outer mgc: 162.4893
## outer mgc: 91.03071
## outer mgc: 98.22766
## outer mgc: 2.172514
## outer mgc: 5.116226
## outer mgc: 6.162891
## outer mgc: 5.261146
## outer mgc: 8.460502
## outer mgc: 3.897932
## outer mgc: 3.234142
## outer mgc: 0.9433855
## outer mgc: 0.3326979
## outer mgc: 0.01294656
## outer mgc: 0.009230682
## outer mgc: 0.001157467
## outer mgc: 0.0001261349
## outer mgc: 1.271994e-05
Extracting estimates and their SEs and showing results
best <- model$env$last.par.best
rep <- sdreport(model)
## outer mgc: 1.271994e-05
## outer mgc: 0.017538
## outer mgc: 0.01756344
## outer mgc: 0.1027843
## outer mgc: 0.1028098
## outer mgc: 0.01999046
## outer mgc: 0.02000952
print(best)
## b0 b1 logSigma
## 4.0781824 1.9090907 0.3451266
print(rep)
## sdreport(.) result
## Estimate Std. Error
## b0 4.0781824 0.7039414
## b1 1.9090907 0.1554747
## logSigma 0.3451266 0.2236068
## Maximum gradient component: 1.271994e-05