A review of statistical computing with TMB

class: center, middle, inverse, title-slide

.title[
# A review of statistical computing with TMB
]
.author[
### Andrea Havron<br>NOAA Fisheries, OST
]

---

layout: true

.footnote[U.S. Department of Commerce | National Oceanic and Atmospheric Administration | National Marine Fisheries Service]

code.cpp{
  font-size: 14px;
}
code.r{
  font-size: 14px;
}

</style>

---
# Template Model Builder
<img src="data:image/png;base64,#static/tmb-wordcloud.png" width="45%" style="display: block; margin: auto;" />

Wordcloud generated from: ([Kristensen et al., 2015](https://backend.orbit.dtu.dk/ws/portalfiles/portal/123599315/Publishers_version.pdf))
---
# Statistical Computing Review

---
# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. **Specify the model** <br><br>
`\(y ~ \sim Binomial(n, p)\)`
]
---
# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. **Calculate the likelihood**<br><br>
`\begin{align}
L(p; n, y) &= \frac{n!}{y!(n-y)!}p^y(1-p)^{n-y}
\end{align}`
]

.pull-right[
`\(y = dbinom(x = 30, n = 100, p)\)`
<div id="htmlwidget-65cc438f8d7e14c2fc7a" style="width:450px;height:375px;" class="plotly html-widget "></div>
<script type="application/json" data-for="htmlwidget-65cc438f8d7e14c2fc7a">{"x":{"visdat":{"516c3c1f6cb1":["function () ","plotlyVisDat"]},"cur_data":"516c3c1f6cb1","attrs":{"516c3c1f6cb1":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"y"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001,0.001,0.003,0.005,0.009,0.015,0.024,0.034,0.046,0.058,0.069,0.079,0.085,0.087,0.085,0.079,0.071,0.06,0.049,0.039,0.029,0.021,0.015,0.01,0.006,0.004,0.002,0.001,0.001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---

# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. **Calculate the negative log-likelihood**

`\begin{align}
-\ell(p; n, y) = &-[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)\\
&+ (n-y)ln(1-p)]
\end{align}`
]

.pull-right[
`\(nll = -log(dbinom(x = 30, n = 100, p))\)`
<div id="htmlwidget-cb1d5d2496ce258c9121" style="width:450px;height:375px;" class="plotly html-widget "></div>
<script type="application/json" data-for="htmlwidget-cb1d5d2496ce258c9121">{"x":{"visdat":{"516cf34313c":["function () ","plotlyVisDat"]},"cur_data":"516cf34313c","attrs":{"516cf34313c":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"nll"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[null,80.217,60.133,48.687,40.782,34.82,30.092,26.216,22.966,20.198,17.811,15.734,13.914,12.313,10.899,9.648,8.54,7.56,6.693,5.93,5.261,4.678,4.174,3.744,3.382,3.084,2.847,2.668,2.542,2.468,2.444,2.468,2.537,2.651,2.808,3.007,3.248,3.528,3.848,4.207,4.604,5.04,5.514,6.025,6.575,7.162,7.787,8.45,9.152,9.893,10.673,11.493,12.354,13.256,14.201,15.189,16.221,17.299,18.425,19.599,20.823,22.099,23.43,24.817,26.262,27.769,29.34,30.979,32.688,34.473,36.336,38.284,40.321,42.453,44.686,47.029,49.489,52.076,54.801,57.675,60.713,63.931,67.347,70.985,74.869,79.032,83.51,88.351,93.611,99.363,105.7,112.743,120.66,129.683,140.153,152.598,167.904,187.731,215.806,264.021,null],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---

# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
`\begin{align}
-\ell(p; n, y) = &-[ln\big(\frac{n!}{y!(n-y)!}\big) + yln(p)\\
&+ (n-y)ln(1-p)]
\end{align}`

4. **Calculate the derivative w.r.t. `\(p\)`**
`\begin{align}
\frac{d(\ell(p; n, y))}{dp} &= \frac{y}{p}- \frac{n-y}{1-p}
\end{align}`
]

.pull-right[
`\(nll = -log(dbinom(x = 30, n = 100, p))\)`
<div id="htmlwidget-4ee9a0d23a5d92cf8569" style="width:450px;height:375px;" class="plotly html-widget "></div>
<script type="application/json" data-for="htmlwidget-4ee9a0d23a5d92cf8569">{"x":{"visdat":{"516c2ef07e23":["function () ","plotlyVisDat"]},"cur_data":"516c2ef07e23","attrs":{"516c2ef07e23":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"nll"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[null,80.217,60.133,48.687,40.782,34.82,30.092,26.216,22.966,20.198,17.811,15.734,13.914,12.313,10.899,9.648,8.54,7.56,6.693,5.93,5.261,4.678,4.174,3.744,3.382,3.084,2.847,2.668,2.542,2.468,2.444,2.468,2.537,2.651,2.808,3.007,3.248,3.528,3.848,4.207,4.604,5.04,5.514,6.025,6.575,7.162,7.787,8.45,9.152,9.893,10.673,11.493,12.354,13.256,14.201,15.189,16.221,17.299,18.425,19.599,20.823,22.099,23.43,24.817,26.262,27.769,29.34,30.979,32.688,34.473,36.336,38.284,40.321,42.453,44.686,47.029,49.489,52.076,54.801,57.675,60.713,63.931,67.347,70.985,74.869,79.032,83.51,88.351,93.611,99.363,105.7,112.743,120.66,129.683,140.153,152.598,167.904,187.731,215.806,264.021,null],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---

# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
4. Calculate the derivate wrt p
5. **Set to `\(0\)` and solve for MLE**<br>

`\begin{align}
0 &= \frac{y}{p}- \frac{n-y}{1-p} \\
E[p] &= \frac{y}{n} \\
E[y] &= np
\end{align}`

]

.pull-right[
`\(nll = -log(dbinom(x = 30, n = 100, p))\)`
<div id="htmlwidget-62db5c7fff9c0f9712d9" style="width:450px;height:375px;" class="plotly html-widget "></div>
<script type="application/json" data-for="htmlwidget-62db5c7fff9c0f9712d9">{"x":{"visdat":{"516c7b127c83":["function () ","plotlyVisDat"]},"cur_data":"516c7b127c83","attrs":{"516c7b127c83":{"x":{},"y":{},"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"},"516c7b127c83.1":{"x":[0.1,0.5],"y":[2.44433456453297,2.44433456453297],"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter","line":{"color":"red"},"showlegend":false,"inherit":true}},"layout":{"width":450,"height":375,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":"p"},"yaxis":{"domain":[0,1],"automargin":true,"title":"nll"},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[null,80.217,60.133,48.687,40.782,34.82,30.092,26.216,22.966,20.198,17.811,15.734,13.914,12.313,10.899,9.648,8.54,7.56,6.693,5.93,5.261,4.678,4.174,3.744,3.382,3.084,2.847,2.668,2.542,2.468,2.444,2.468,2.537,2.651,2.808,3.007,3.248,3.528,3.848,4.207,4.604,5.04,5.514,6.025,6.575,7.162,7.787,8.45,9.152,9.893,10.673,11.493,12.354,13.256,14.201,15.189,16.221,17.299,18.425,19.599,20.823,22.099,23.43,24.817,26.262,27.769,29.34,30.979,32.688,34.473,36.336,38.284,40.321,42.453,44.686,47.029,49.489,52.076,54.801,57.675,60.713,63.931,67.347,70.985,74.869,79.032,83.51,88.351,93.611,99.363,105.7,112.743,120.66,129.683,140.153,152.598,167.904,187.731,215.806,264.021,null],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null},{"x":[0.1,0.5],"y":[2.44433456453297,2.44433456453297],"mode":"lines","type":"scatter","line":{"color":"red"},"showlegend":false,"marker":{"color":"rgba(255,127,14,1)","line":{"color":"rgba(255,127,14,1)"}},"error_y":{"color":"rgba(255,127,14,1)"},"error_x":{"color":"rgba(255,127,14,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]

---
# Maximum Likelihood Inference

#### What is the likelihood of getting 30 heads in 100 coin flips?

.pull-left[
1. Specify the model
2. Calculate the likelihood
3. Calculate the negative log-likelihood
4. Calculate the derivate wrt p
5. Set to `\(0\)` and solve for MLE
`\begin{align}
0 &= \frac{y}{p}- \frac{n-y}{1-p}, \hspace{2mm} E[y] = np \\
\end{align}`
6. **Approximate `\(Var[p]\)` using the second derivative**<br>
`\begin{align}
&-\frac{y}{p^2} - \frac{(n-y)}{(1-p)^2} = -\frac{np}{p^2} - \frac{(n-np)}{(1-p)^2}\\
 &= -\frac{n}{p} - \frac{n}{1-p}\\
l''(p) &= -\frac{n(1-p)}{p}
\end{align}`
]

.pull-right[
`\begin{align}
Var[p] &\approx -\frac{1}{l''(p)} \approx \frac{p(1-p)}{n}\\
SE[p] &\approx \sqrt{ .3(.7)/100} \approx .00458
\end{align}`
<div id="htmlwidget-9daccd87cff7a4d88ed6" style="width:360px;height:300px;" class="plotly html-widget "></div>
<script type="application/json" data-for="htmlwidget-9daccd87cff7a4d88ed6">{"x":{"visdat":{"516c7ecb434e":["function () ","plotlyVisDat"]},"cur_data":"516c7ecb434e","attrs":{"516c7ecb434e":{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001,0.001,0.003,0.005,0.009,0.015,0.024,0.034,0.046,0.058,0.069,0.079,0.085,0.087,0.085,0.079,0.071,0.06,0.049,0.039,0.029,0.021,0.015,0.01,0.006,0.004,0.002,0.001,0.001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter"},"516c7ecb434e.1":{"x":0.3,"y":0.0867838647534281,"mode":"markers","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter","showlegend":false,"inherit":true},"516c7ecb434e.2":{"x":0.208348486100883,"y":0.04,"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter","xend":0.208348486100883,"yend":0,"line":{"color":"red"},"showlegend":false,"inherit":true},"516c7ecb434e.3":{"x":0.391651513899117,"y":0.04,"mode":"lines","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter","xend":0.391651513899117,"yend":0,"line":{"color":"red"},"showlegend":false,"inherit":true},"516c7ecb434e.4":{"x":0.31,"y":0.01,"mode":"text","alpha_stroke":1,"sizes":[10,100],"spans":[1,20],"type":"scatter","text":"95% <br />CI","showlegend":false,"inherit":true}},"layout":{"width":360,"height":300,"margin":{"b":40,"l":60,"t":25,"r":10},"xaxis":{"domain":[0,1],"automargin":true,"title":[]},"yaxis":{"domain":[0,1],"automargin":true,"title":[]},"hovermode":"closest","showlegend":false},"source":"A","config":{"modeBarButtonsToAdd":["hoverclosest","hovercompare"],"showSendToCloud":false},"data":[{"x":[0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09,0.1,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18,0.19,0.2,0.21,0.22,0.23,0.24,0.25,0.26,0.27,0.28,0.29,0.3,0.31,0.32,0.33,0.34,0.35,0.36,0.37,0.38,0.39,0.4,0.41,0.42,0.43,0.44,0.45,0.46,0.47,0.48,0.49,0.5,0.51,0.52,0.53,0.54,0.55,0.56,0.57,0.58,0.59,0.6,0.61,0.62,0.63,0.64,0.65,0.66,0.67,0.68,0.69,0.7,0.71,0.72,0.73,0.74,0.75,0.76,0.77,0.78,0.79,0.8,0.81,0.82,0.83,0.84,0.85,0.86,0.87,0.88,0.89,0.9,0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1],"y":[0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0.001,0.001,0.003,0.005,0.009,0.015,0.024,0.034,0.046,0.058,0.069,0.079,0.085,0.087,0.085,0.079,0.071,0.06,0.049,0.039,0.029,0.021,0.015,0.01,0.006,0.004,0.002,0.001,0.001,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0],"mode":"lines","type":"scatter","marker":{"color":"rgba(31,119,180,1)","line":{"color":"rgba(31,119,180,1)"}},"error_y":{"color":"rgba(31,119,180,1)"},"error_x":{"color":"rgba(31,119,180,1)"},"line":{"color":"rgba(31,119,180,1)"},"xaxis":"x","yaxis":"y","frame":null},{"x":[0.3],"y":[0.0867838647534281],"mode":"markers","type":"scatter","showlegend":false,"marker":{"color":"rgba(255,127,14,1)","line":{"color":"rgba(255,127,14,1)"}},"error_y":{"color":"rgba(255,127,14,1)"},"error_x":{"color":"rgba(255,127,14,1)"},"line":{"color":"rgba(255,127,14,1)"},"xaxis":"x","yaxis":"y","frame":null},{"x":[0.208348486100883,0.208348486100883],"y":[0.04,0],"mode":"lines","type":"scatter","line":{"color":"red"},"showlegend":false,"marker":{"color":"rgba(44,160,44,1)","line":{"color":"rgba(44,160,44,1)"}},"error_y":{"color":"rgba(44,160,44,1)"},"error_x":{"color":"rgba(44,160,44,1)"},"xaxis":"x","yaxis":"y","frame":null},{"x":[0.391651513899117,0.391651513899117],"y":[0.04,0],"mode":"lines","type":"scatter","line":{"color":"red"},"showlegend":false,"marker":{"color":"rgba(214,39,40,1)","line":{"color":"rgba(214,39,40,1)"}},"error_y":{"color":"rgba(214,39,40,1)"},"error_x":{"color":"rgba(214,39,40,1)"},"xaxis":"x","yaxis":"y","frame":null},{"x":[0.31],"y":[0.01],"mode":"text","type":"scatter","text":"95% <br />CI","showlegend":false,"marker":{"color":"rgba(148,103,189,1)","line":{"color":"rgba(148,103,189,1)"}},"error_y":{"color":"rgba(148,103,189,1)"},"error_x":{"color":"rgba(148,103,189,1)"},"line":{"color":"rgba(148,103,189,1)"},"xaxis":"x","yaxis":"y","frame":null}],"highlight":{"on":"plotly_click","persistent":false,"dynamic":false,"selectize":false,"opacityDim":0.2,"selected":{"opacity":1},"debounce":0},"shinyEvents":["plotly_hover","plotly_click","plotly_selected","plotly_relayout","plotly_brushed","plotly_brushing","plotly_clickannotation","plotly_doubleclick","plotly_deselect","plotly_afterplot","plotly_sunburstclick"],"base_url":"https://plot.ly"},"evals":[],"jsHooks":[]}</script>
]
---
# Multivariate asymptotics

* The second derivative measures the curvature of the likelihood surface
* For models with N parameters, the curvature is represented by an NxN **Hessian** matrix of 2nd partial derivatives
* Inverting the negative Hessian gives us a covariance matrix

.pull-left[
`\begin{align}
(\mathbb{H}_{f})_{i,j} &= \frac{\partial^2f}{\partial \theta_{i}, \partial x\theta_{j}} = \frac{-1}{Var(\Theta)}
\end{align}`
<br>
[**What causes a singular covariance matrix?**](https://andrea-havron.shinyapps.io/mvnplot/)
]

.pull-right[
<img src="data:image/png;base64,#00_00_Statistical_Computing_Review_files/figure-html/mvnorm-1.png" width="80%" />
]

---
# Automatic Differentiation

<img src="data:image/png;base64,#static/AD_triangle.png" width="65%" style="display: block; margin: auto;" />
---
# What is AD?
<br>
<table>
 <thead>
  <tr>
   <th style="text-align:left;">  </th>
   <th style="text-align:left;"> Automatic Differentiation </th>
   <th style="text-align:left;"> Symbolic Differentiation </th>
   <th style="text-align:left;"> Numerical Differentiation </th>
  </tr>
 </thead>
<tbody>
  <tr>
   <td style="text-align:left;">  </td>
   <td style="text-align:left;"> Derivatives calculated automatically using the chain rule </td>
   <td style="text-align:left;"> Computer program converts function into exact derivative function </td>
   <td style="text-align:left;"> Approximation that relies on finite differences </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Efficiency: </td>
   <td style="text-align:left;"> forward mode: O(n); reverse mode: O(m) </td>
   <td style="text-align:left;"> O(2n) </td>
   <td style="text-align:left;"> O(n) </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Accuracy: </td>
   <td style="text-align:left;"> Machine Precision </td>
   <td style="text-align:left;"> Exact </td>
   <td style="text-align:left;"> Trade-off between truncation error and round-off error </td>
  </tr>
  <tr>
   <td style="text-align:left;"> Higher order derivatives: </td>
   <td style="text-align:left;"> Easy </td>
   <td style="text-align:left;"> Difficult due to complexity </td>
   <td style="text-align:left;"> Difficult due to error accumulation </td>
  </tr>
</tbody>
</table>
*O(n): Order of operations is based on the number of parameters <br>
*O(m): Order of operations is based on the number of outputs, in MLE inference, this is the nll
---
# Computational Graph (Tape)
<br>
.three-column-left[

```cpp
//Program
v1: x = ?
v2: y = ?
v3: a = x * y
v4: b = sin(y)
v5: z = a + b
```
]
.three-column[
<img src="data:image/png;base64,#static/comp-graph.png" width="100%" style="display: block; margin: auto auto auto 0;" />
]
.three-column-left[

```cpp
//Reverse Mode
dz = ?
da = dz
db = dz
dx = yda
dy = xda + cos(y)db
```
]

---
# Reverse Mode
.pull-left[
**Static (TMB: CppAD, TMBad)**<br>
The graph is constructed once before execution
.p[
- Less flexibility with conditional statements that depend on parameters. 
- Atomic functions can be used when conditional statements depend on parameters
- High portability 
- Graph optimization possible
]]

.pull-right[
**Dynamic (Stan: Stan Math Library, ADMB: AUTODIF)**<br>
The graph is defined as forward computation is executed at every iteration
.p[
- Flexibility with conditional statements
- Optimization routine implemented into executable
- Less time for graph optimization
]
]

---

# TMB AD Systems
<br>
.pull-left[
**CppAD**
- [CppAD package](https://coin-or.github.io/CppAD/doc/cppad.htm)
- Developed and maintained by the CppAD community

```r
TMB::compile("mymodel.cpp")
TMB::compile("mymodel.cpp", framework = "CppAD")
```

]
.pull-right[
**TMBad**<br>
- TMBad is available with TMB 1.8.0 and higher
- Improved efficiencies over CppAD
- Developed and maintained by Kasper Kristensen

```r
TMB::compile("mymodel.cpp", framework = "TMBad")
```

]
---
# R and C++ Type Systems

<img src="data:image/png;base64,#static/Type_Triangle.png" width="65%" style="display: block; margin: auto;" />
---
# Dynamic vs. Static Typing
<br>
.bluebox[
**Type System** <br>
A logical system comprising a set of rules that assigns a property called a type (for example, integer, floating point, string) to every "term" (a word, phrase, or other set of symbols) <br>
Source: [Wikipedia](https://en.wikipedia.org/wiki/Type_system)
]

.pull-left[
**R: Dynamic**
.p[
- Type checking occurs at run time
- The values and types associated with names can change
- Change in type tends to be implicit
]]

.pull-right[
**C++: Static**
.p[
- Type checking occurs at compile time
- The values associated with a given name can be limited to just a few types and may be unchangeable
- Change in type tends to be explicit
]]
---
# Dynamic vs. Static Typing
<br>
.pull-left[
**R: Dynamic**
- Flexible 
- Easy to debug at runtime
- Code is easier to read and write
- More errors detected later in development
- More tests are needed to detect errors
- Slower to run

```r
`+`  <- `-`
1 + 1
```

```
## [1] 0
```

]

.pull-right[
**C++: Static**
- Errors detected early in development
- Fewer errors occur at runtime
- Optimized code runs faster
- Code is harder to read and write
- Difficult to debug
<br><br>
<img src="data:image/png;base64,#static/RStudioBomb.png" width="65%" style="display: block; margin: auto;" />

]

---
# What is the Template in TMB?
<br>
**Templated C++**
<br>
* Generic programming
* Allows developer to write functions and classes that are independent of Type
* Templates are expanded at compile time

.pull-left[

```cpp

template <class Type>
  Type add(Type x, Type y){
  return x + y;
}

int main(){
  int a = 1;
  int b = 2;
  double c = 1.1;
  double d = 2.1;
  int d = add(a,b);
  double e = add(c,d);
}

```
]
.pull-right[

```cpp
int add(int x, int y){
  return x + y;
}
double add(double x, double y){
  return x + y;
}
```
]