<!DOCTYPE html>
<!-- Generated by pkgdown: do not edit by hand --><html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=UTF-8">
<meta charset="utf-8">
<meta http-equiv="X-UA-Compatible" content="IE=edge">
<meta name="viewport" content="width=device-width, initial-scale=1.0">
<title>hmclearn package:  Linear Regression Example • hmclearn</title>
<!-- jquery --><script src="https://cdnjs.cloudflare.com/ajax/libs/jquery/3.4.1/jquery.min.js" integrity="sha256-CSXorXvZcTkaix6Yvo6HppcZGetbYMGWSFlBw8HfCJo=" crossorigin="anonymous"></script><!-- Bootstrap --><link href="https://cdnjs.cloudflare.com/ajax/libs/bootswatch/3.4.0/sandstone/bootstrap.min.css" rel="stylesheet" crossorigin="anonymous">
<script src="https://cdnjs.cloudflare.com/ajax/libs/twitter-bootstrap/3.4.1/js/bootstrap.min.js" integrity="sha256-nuL8/2cJ5NDSSwnKD8VqreErSWHtnEP9E7AySL+1ev4=" crossorigin="anonymous"></script><!-- bootstrap-toc --><link rel="stylesheet" href="../bootstrap-toc.css">
<script src="../bootstrap-toc.js"></script><!-- Font Awesome icons --><link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/all.min.css" integrity="sha256-mmgLkCYLUQbXn0B1SRqzHar6dCnv9oZFPEC1g1cwlkk=" crossorigin="anonymous">
<link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/font-awesome/5.12.1/css/v4-shims.min.css" integrity="sha256-wZjR52fzng1pJHwx4aV2AO3yyTOXrcDW7jBpJtTwVxw=" crossorigin="anonymous">
<!-- clipboard.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/clipboard.js/2.0.6/clipboard.min.js" integrity="sha256-inc5kl9MA1hkeYUt+EC3BhlIgyp/2jDIyBLS6k3UxPI=" crossorigin="anonymous"></script><!-- headroom.js --><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/headroom.min.js" integrity="sha256-AsUX4SJE1+yuDu5+mAVzJbuYNPHj/WroHuZ8Ir/CkE0=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/headroom/0.11.0/jQuery.headroom.min.js" integrity="sha256-ZX/yNShbjqsohH1k95liqY9Gd8uOiE1S4vZc+9KQ1K4=" crossorigin="anonymous"></script><!-- pkgdown --><link href="../pkgdown.css" rel="stylesheet">
<script src="../pkgdown.js"></script><meta property="og:title" content="hmclearn package:  Linear Regression Example">
<meta property="og:description" content="hmclearn">
<!-- mathjax --><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js" integrity="sha256-nvJJv9wWKEm88qvoQl9ekL2J+k/RWIsaSScxxlsrv8k=" crossorigin="anonymous"></script><script src="https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/config/TeX-AMS-MML_HTMLorMML.js" integrity="sha256-84DKXVJXs0/F8OTMzX4UR909+jtl4G7SPypPavF+GfA=" crossorigin="anonymous"></script><!--[if lt IE 9]>
<script src="https://oss.maxcdn.com/html5shiv/3.7.3/html5shiv.min.js"></script>
<script src="https://oss.maxcdn.com/respond/1.4.2/respond.min.js"></script>
<![endif]-->
</head>
<body data-spy="scroll" data-target="#toc">
    <div class="container template-article">
      <header><div class="navbar navbar-default navbar-fixed-top" role="navigation">
  <div class="container">
    <div class="navbar-header">
      <button type="button" class="navbar-toggle collapsed" data-toggle="collapse" data-target="#navbar" aria-expanded="false">
        <span class="sr-only">Toggle navigation</span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
        <span class="icon-bar"></span>
      </button>
      <span class="navbar-brand">
        <a class="navbar-link" href="../index.html">hmclearn</a>
        <span class="version label label-default" data-toggle="tooltip" data-placement="bottom" title="Released version">0.0.2</span>
      </span>
    </div>

    <div id="navbar" class="navbar-collapse collapse">
      <ul class="nav navbar-nav">
<li>
  <a href="../index.html">
    <span class="fas fa fas fa-home fa-lg"></span>
     
  </a>
</li>
<li>
  <a href="../reference/index.html">Reference</a>
</li>
<li class="dropdown">
  <a href="#" class="dropdown-toggle" data-toggle="dropdown" role="button" aria-expanded="false">
    Articles
     
    <span class="caret"></span>
  </a>
  <ul class="dropdown-menu" role="menu">
<li>
      <a href="../articles/linear_mixed_effects_hmclearn.html">hmclearn package:  Linear Mixed Effects Regression Example</a>
    </li>
    <li>
      <a href="../articles/linear_regression_hmclearn.html">hmclearn package:  Linear Regression Example</a>
    </li>
    <li>
      <a href="../articles/logistic_mixed_effects_hmclearn.html">hmclearn package:  Logistic Mixed Effects Regression Example</a>
    </li>
    <li>
      <a href="../articles/logistic_regression_hmclearn.html">hmclearn:  Logistic Regression Example</a>
    </li>
    <li>
      <a href="../articles/poisson_regression_hmclearn.html">hmclearn:  Poisson Regression Example</a>
    </li>
  </ul>
</li>
<li>
  <a href="../news/index.html">Changelog</a>
</li>
      </ul>
<ul class="nav navbar-nav navbar-right"></ul>
</div>
<!--/.nav-collapse -->
  </div>
<!--/.container -->
</div>
<!--/.navbar -->

      

      </header><div class="row">
  <div class="col-md-9 contents">
    <div class="page-header toc-ignore">
      <h1 data-toc-skip>hmclearn package: Linear Regression Example</h1>
                        <h4 class="author">Samuel Thomas</h4>
            
            <h4 class="date">2020-06-01</h4>
      
      
      <div class="hidden name"><code>linear_regression_hmclearn.Rmd</code></div>

    </div>

    
    
<div id="introduction" class="section level2">
<h2 class="hasAnchor">
<a href="#introduction" class="anchor"></a>Introduction</h2>
<p>This vignette demonstrates fitting a linear regression model via Hamiltonian Monte Carlo (HMC) using the <strong>hmclearn</strong> package.</p>
<p><span class="math display">\[
\begin{aligned}
\mathbf{y} &amp;= \mathbf{X}\boldsymbol\beta + \boldsymbol\epsilon \\
\boldsymbol\epsilon &amp;\sim N(0, \sigma_\epsilon^2)
\end{aligned}
\]</span></p>
<p>HMC requires the specification of the log posterior to a proportional constant. In addition, HMC uses the gradient of the log posterior to guide simulations. The full derivations are provided here.</p>
<p>The vector of responses is <span class="math inline">\(\mathbf{y} = (y_1, ..., y_n)^T\)</span>. The covariate values for the th subject are <span class="math inline">\(\mathbf{x}_i^T = (x_{i0}, ..., x_{iq})\)</span> for <span class="math inline">\(q\)</span> covariates plus an intercept. We write the full design matrix as <span class="math inline">\(\mathbf{X} = (\mathbf{x}_1^T, ..., \mathbf{x}_n^T)^T \in \mathbb{R}^{n\times(q+1)}\)</span> for <span class="math inline">\(n\)</span> observations. The regression coefficients are a vector of length <span class="math inline">\(q + 1\)</span>, <span class="math inline">\(\boldsymbol\beta = (\beta_0, ..., \beta_q)^T\)</span>. The error term for each subject is <span class="math inline">\(\epsilon_i\)</span>. All error terms <span class="math inline">\(\boldsymbol\epsilon = (\epsilon_1, ..., \epsilon_n)^T\)</span> are assumed to be independent and normally distributed with mean zero and constant variance <span class="math inline">\(\sigma_\epsilon^2\)</span>.</p>
</div>
<div id="derive-log-posterior-and-gradient-for-hmc" class="section level2">
<h2 class="hasAnchor">
<a href="#derive-log-posterior-and-gradient-for-hmc" class="anchor"></a>Derive log posterior and gradient for HMC</h2>
<p>We specify the likelihood function for linear regression</p>
<p><span class="math display">\[
f(\mathbf{y} | X, \boldsymbol\beta, \sigma_\epsilon^2) = \frac{1}{(2\pi\sigma_\epsilon^2)^{n/2}}\exp{\left(-\frac{1}{2\sigma_\epsilon^2} (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) \right)}.
\]</span></p>
<p>The posterior is defined with priors for <span class="math inline">\(\boldsymbol\beta\)</span> and <span class="math inline">\(\sigma_\epsilon^2\)</span>, such that</p>
<p><span class="math display">\[
\begin{aligned}
f(\boldsymbol\beta, \sigma_\epsilon^2| \mathbf{y}, \mathbf{X}) &amp;\propto f(\mathbf{y} |\mathbf{X}, \boldsymbol\beta, \sigma_\epsilon^2) \pi(\boldsymbol\beta, \sigma_\epsilon^2), \\
&amp;\propto f(\mathbf{y} | \mathbf{X}, \boldsymbol\beta, \sigma_\epsilon^2) \pi(\boldsymbol\beta) \pi(\sigma_\epsilon^2).
\end{aligned}
\]</span></p>
<p>We specify a multivariate normal prior for <span class="math inline">\(\boldsymbol\beta\)</span> and Inverse Gamma for <span class="math inline">\(\sigma_\epsilon^2\)</span>. Note that Inverse Gamma can be problematic in certain models where the variance is low (Gelman 2006). The prior distributions in this example are</p>
<p><span class="math display">\[
\begin{aligned}
\pi(\boldsymbol\beta | \sigma_\beta^2) &amp;\propto N(0, \sigma_\beta^2 \mathbf{I}), \\
\pi(\sigma_\epsilon^2 | a, b) &amp;\sim IG(a, b).
\end{aligned}
\]</span></p>
<p>The pdf for <span class="math inline">\(\sigma_\epsilon^2\)</span> prior is</p>
<p><span class="math display">\[
f(\sigma_\epsilon^2) = \frac{b^a}{\Gamma(a)} (\sigma^2)^{-a-1} \exp\left(-\frac{b}{\sigma_\epsilon^2} \right).
\]</span></p>
<p>We perform a log-transformation of <span class="math inline">\(\sigma^2\)</span> to allow proposals across entire real number line. For unconstrained HMC simulations, all parameters must support all real numbers. We provide details for the transformation of variable,</p>
<p><span class="math display">\[
\begin{aligned}
\gamma &amp;= \log \sigma_\epsilon^2, \\
\sigma_\epsilon^2 &amp;= g^{-1}(\gamma) = e^\gamma, \\
\pi_\gamma(\gamma | a, b) &amp;= p_{\sigma_\epsilon^2}(g^{-1}(\gamma))\left\lvert \frac{d\sigma_\epsilon^2}{d\gamma} \right\rvert, \\
&amp;= \frac{b^a}{\Gamma(a)}(e^\gamma)^{-a-1} \exp\left(-\frac{b}{\sigma^2} \right) \left\lvert e^\gamma  \right\rvert, \\
&amp;= \frac{b^a}{\Gamma(a)}e^{-a\gamma}\exp\left(-\frac{b}{e^\gamma} \right) \\
&amp;\propto e^{-a\gamma}\exp\left(-\frac{b}{e
^\gamma} \right), \\
\log \pi(\gamma | a, b) &amp;\propto -a\gamma - be^{-\gamma}.
\end{aligned}
\]</span></p>
<p>The posterior and log posterior for linear regression, based on <span class="math inline">\(\beta\)</span> and <span class="math inline">\(\gamma\)</span>, can now be derived.</p>
<p><span class="math display">\[
\begin{aligned}
f(\boldsymbol\beta, \gamma| \mathbf{y}, \mathbf{X}, \sigma_\beta^2, a, b) &amp;\propto f(\mathbf{y} | \mathbf{X}, \boldsymbol\beta, \gamma) \pi(\boldsymbol\beta, \gamma | \sigma_\beta^2, a, b), \\
&amp;\propto f(\mathbf{y}| \mathbf{X}, \boldsymbol\beta, \gamma) \pi(\boldsymbol\beta | \sigma_\beta^2) \pi(\gamma | a, b), \\
&amp;\propto \frac{e^{-\gamma n/2}}{(2\pi)^{n/2}}\exp{\left(-\frac{e^{-\gamma}}{2} (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) \right)} e^{-\frac{1}{2}\boldsymbol\beta^T\boldsymbol\beta / \sigma_\beta^2} e^{-a\gamma} e^{-b e^{-\gamma}},  \\
\log f(\boldsymbol\beta, \gamma | \mathbf{y}, \mathbf{X}, \sigma_\beta^2, a, b) &amp;\propto -\frac{\gamma n}{2} -\frac{e^{-\gamma}}{2} (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) - a\gamma - b e^{-\gamma}- \frac{\boldsymbol\beta^T \boldsymbol\beta}{2\sigma_\beta^2},   \\
&amp;\propto -\left(\frac{n}{2} + a \right)\gamma  -\frac{e^{-\gamma}}{2} (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta)- b e^{-\gamma} - \frac{\boldsymbol\beta^T \boldsymbol\beta}{2\sigma_\beta^2}.
\end{aligned}
\]</span></p>
<p>Next, we derive the gradient of the log posterior</p>
<p><span class="math display">\[
\begin{aligned}
\log f(\boldsymbol\beta, \gamma | \mathbf{y}, \mathbf{X}, \sigma_\beta^2, a, b) &amp;\propto -\frac{\gamma n}{2} -\frac{e^{-\gamma}}{2} (\mathbf{y} - \mathbf{X}\boldsymbol\beta)^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) - a\gamma - b e^{-\gamma} - \frac{\boldsymbol\beta^T \boldsymbol\beta}{2\sigma_\beta^2},  \\
&amp;\propto -\left(\frac{n}{2} + a \right)\gamma  -\frac{e^{-\gamma}}{2} (y - X\beta)^T(y-X\beta)- b e^{-\gamma}- \frac{\boldsymbol\beta^T \boldsymbol\beta}{2\sigma_\beta^2},  \\
\nabla_{\boldsymbol\beta }\log f(\boldsymbol\beta, \gamma | \mathbf{y}, \mathbf{X}, \sigma_\beta^2, a, b) &amp;\propto -\frac{e^{-\gamma}}{2}(-2)\mathbf{X}^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) - \boldsymbol\beta/\sigma_\beta^2, \\
&amp;\propto e^{-\gamma} \mathbf{X}^T(\mathbf{y}-\mathbf{X}\boldsymbol\beta) - \boldsymbol\beta/\sigma_\beta^2, \\
\nabla_\gamma\log f(\beta, \gamma | y, \mathbf{X}, \sigma_\beta^2, a, b) &amp;\propto -\left(\frac{n}{2} + a \right) + \frac{e^{-\gamma}}{2} (y - X\beta)^T(y-X\beta)+ b e^{-\gamma}. 
\end{aligned}
\]</span></p>
</div>
<div id="linear-regression-example-data" class="section level2">
<h2 class="hasAnchor">
<a href="#linear-regression-example-data" class="anchor"></a>Linear Regression Example Data</h2>
<div class="sourceCode" id="cb1"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/library.html">library</a></span>(<span class="no">hmclearn</span>)</pre></body></html></div>
<p>The user must define provide the design matrix directly for use in <strong>hmclearn</strong>. Our first step is to load the data and store the design matrix <span class="math inline">\(X\)</span> and dependent variable vector <span class="math inline">\(y\)</span>.</p>
<p>Load Scots Race data (Atkinson 1986, Venables and Ripley 2002) and display the first few rows of the design matrix <span class="math inline">\(X\)</span>. This example also appears in Agresti (2015), and we compare results to his.</p>
<div class="sourceCode" id="cb2"><html><body><pre class="r"><span class="no">ScotsRaces</span> <span class="kw">&lt;-</span> <span class="kw pkg">MASS</span><span class="kw ns">::</span><span class="no"><a href="https://rdrr.io/pkg/MASS/man/hills.html">hills</a></span>

<span class="no">X</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span>(<span class="fl">1</span>, <span class="fu"><a href="https://rdrr.io/r/base/matrix.html">as.matrix</a></span>(<span class="no">ScotsRaces</span>[, -<span class="fu"><a href="https://rdrr.io/r/base/which.html">which</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">ScotsRaces</span>) <span class="kw">==</span> <span class="st">"time"</span>)]))

<span class="co"># add interaction</span>
<span class="no">X</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/cbind.html">cbind</a></span>(<span class="no">X</span>, <span class="no">X</span>[, <span class="st">"dist"</span>] * <span class="no">X</span>[, <span class="st">"climb"</span>])
<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>)[<span class="fu"><a href="https://rdrr.io/r/base/nrow.html">ncol</a></span>(<span class="no">X</span>)] <span class="kw">&lt;-</span> <span class="st">"climb_distance"</span>
<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>)[<span class="fl">1</span>] <span class="kw">&lt;-</span> <span class="st">"intercept"</span>

<span class="no">y</span> <span class="kw">&lt;-</span> <span class="no">ScotsRaces</span>$<span class="no">time</span>

<span class="fu"><a href="https://rdrr.io/r/utils/head.html">head</a></span>(<span class="no">X</span>, <span class="fl">3</span>)
<span class="co">#&gt;              intercept dist climb climb_distance</span>
<span class="co">#&gt; Greenmantle          1  2.5   650           1625</span>
<span class="co">#&gt; Carnethy             1  6.0  2500          15000</span>
<span class="co">#&gt; Craig Dunain         1  6.0   900           5400</span></pre></body></html></div>
</div>
<div id="comparison-model---frequentist" class="section level2">
<h2 class="hasAnchor">
<a href="#comparison-model---frequentist" class="anchor"></a>Comparison model - Frequentist</h2>
<p>To compare results, we first fit a standard linear model using the frequentist function <em>lm</em>. This formula for the linear model omits the automatic inclusion of the intercept. We code the formula in this way since our design matrix already includes the intercept.</p>
<div class="sourceCode" id="cb3"><html><body><pre class="r"><span class="no">f</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/stats/lm.html">lm</a></span>(<span class="no">y</span> ~ <span class="no">X</span>-<span class="fl">1</span>)
<span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">f</span>)
<span class="co">#&gt; </span>
<span class="co">#&gt; Call:</span>
<span class="co">#&gt; lm(formula = y ~ X - 1)</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Residuals:</span>
<span class="co">#&gt;     Min      1Q  Median      3Q     Max </span>
<span class="co">#&gt; -25.994  -4.968  -2.220   2.381  56.115 </span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Coefficients:</span>
<span class="co">#&gt;                   Estimate Std. Error t value Pr(&gt;|t|)    </span>
<span class="co">#&gt; Xintercept       9.3954374  6.8790233   1.366  0.18183    </span>
<span class="co">#&gt; Xdist            4.1489201  0.8352489   4.967 2.36e-05 ***</span>
<span class="co">#&gt; Xclimb          -0.0009710  0.0041648  -0.233  0.81718    </span>
<span class="co">#&gt; Xclimb_distance  0.0009831  0.0003070   3.203  0.00314 ** </span>
<span class="co">#&gt; ---</span>
<span class="co">#&gt; Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1</span>
<span class="co">#&gt; </span>
<span class="co">#&gt; Residual standard error: 12.92 on 31 degrees of freedom</span>
<span class="co">#&gt; Multiple R-squared:  0.9744, Adjusted R-squared:  0.9711 </span>
<span class="co">#&gt; F-statistic: 295.1 on 4 and 31 DF,  p-value: &lt; 2.2e-16</span></pre></body></html></div>
<p>Next, we store the parameter values from the Frequentist fit.</p>
<div class="sourceCode" id="cb4"><html><body><pre class="r"><span class="no">beta.freq</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/numeric.html">as.numeric</a></span>(<span class="fu"><a href="https://rdrr.io/r/stats/coef.html">coefficients</a></span>(<span class="no">f</span>))
<span class="no">sigma2.freq</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/stats/sigma.html">sigma</a></span>(<span class="no">f</span>)^<span class="fl">2</span>
<span class="no">theta.freq</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="no">beta.freq</span>, <span class="fu"><a href="https://rdrr.io/r/base/Log.html">log</a></span>(<span class="no">sigma2.freq</span>))</pre></body></html></div>
<p>The functions for the log posterior and its gradient are coded in the <em>linear_posterior</em> and <em>g_linear_posterior</em> functions, respectively, in <strong>hmclearn</strong>. We use the default hyperparameters in these functions for <span class="math inline">\(\sigma_\beta^2\)</span>, <span class="math inline">\(a\)</span>, and <span class="math inline">\(b\)</span>.</p>
</div>
<div id="fit-model-using-hmc" class="section level2">
<h2 class="hasAnchor">
<a href="#fit-model-using-hmc" class="anchor"></a>Fit model using <em>hmc</em>
</h2>
<p>Next, we fit the linear regression model using HMC. A vector of <em>tuning parameter</em> <span class="math inline">\(\epsilon\)</span> values are specified to align with the data.</p>
<div class="sourceCode" id="cb5"><html><body><pre class="r"><span class="no">N</span> <span class="kw">&lt;-</span> <span class="fl">1e4</span>

<span class="no">eps_vals</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fl">1e-1</span>, <span class="fl">1e-2</span>, <span class="fl">1e-4</span>, <span class="fl">5e-6</span>, <span class="fl">3e-3</span>)

<span class="fu"><a href="https://rdrr.io/r/base/Random.html">set.seed</a></span>(<span class="fl">412</span>)
<span class="no">t1.hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()
 <span class="no">f_hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="../reference/hmc.html">hmc</a></span>(<span class="kw">N</span><span class="kw">=</span><span class="no">N</span>, <span class="kw">theta.init</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/rep.html">rep</a></span>(<span class="fl">0</span>, <span class="fl">5</span>),
                    <span class="kw">epsilon</span> <span class="kw">=</span> <span class="no">eps_vals</span>, <span class="kw">L</span> <span class="kw">=</span> <span class="fl">20</span>,
                    <span class="kw">logPOSTERIOR</span> <span class="kw">=</span> <span class="no">linear_posterior</span>,
                    <span class="kw">glogPOSTERIOR</span> <span class="kw">=</span> <span class="no">g_linear_posterior</span>,
                    <span class="kw">param</span><span class="kw">=</span><span class="fu"><a href="https://rdrr.io/r/base/list.html">list</a></span>(<span class="kw">y</span> <span class="kw">=</span> <span class="no">y</span>, <span class="kw">X</span><span class="kw">=</span><span class="no">X</span>),
                    <span class="kw">varnames</span> <span class="kw">=</span> <span class="fu"><a href="https://rdrr.io/r/base/c.html">c</a></span>(<span class="fu"><a href="https://rdrr.io/r/base/colnames.html">colnames</a></span>(<span class="no">X</span>), <span class="st">"log_sigma2"</span>),
                    <span class="kw">parallel</span><span class="kw">=</span><span class="fl">TRUE</span>, <span class="kw">chains</span><span class="kw">=</span><span class="fl">2</span>)
<span class="no">t2.hmc</span> <span class="kw">&lt;-</span> <span class="fu"><a href="https://rdrr.io/r/base/Sys.time.html">Sys.time</a></span>()
<span class="no">t2.hmc</span> - <span class="no">t1.hmc</span>
<span class="co">#&gt; Time difference of 16.2068 secs</span></pre></body></html></div>
</div>
<div id="mcmc-summary-and-diagnostics" class="section level2">
<h2 class="hasAnchor">
<a href="#mcmc-summary-and-diagnostics" class="anchor"></a>MCMC summary and diagnostics</h2>
<p>The acceptance rate for each of the HMC chains is high for this simple example. More complex applications will often need lower acceptance rates (e.g. between 0.6 and 0.9) for optimal computational efficiency, particularly without an automated tuning algorithm such as NUTS (Hoffman, Gelman 2014).</p>
<div class="sourceCode" id="cb6"><html><body><pre class="r"><span class="no">f_hmc</span>$<span class="no">accept</span>/<span class="no">N</span>
<span class="co">#&gt; [1] 0.9740 0.9761</span></pre></body></html></div>
<p>The posterior quantiles are summarized after removing an initial <em>burnin</em> period. The <span class="math inline">\(\hat{R}\)</span> statistics are close to one, indicating that both HMC chains converged to the same distribution. The <span class="math inline">\(\hat{R}\)</span> statistics provide an indication of convergence. Values close to one indicate that the multiple MCMC chains coverged to the same distribution, while values above 1.1 indicate possible convergence problems. All <span class="math inline">\(\hat{R}\)</span> values in our example are close to one.</p>
<div class="sourceCode" id="cb7"><html><body><pre class="r"><span class="fu"><a href="https://rdrr.io/r/base/summary.html">summary</a></span>(<span class="no">f_hmc</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">2000</span>)
<span class="co">#&gt; Summary of MCMC simulation</span>
<span class="co">#&gt;                         2.5%            5%           25%          50%</span>
<span class="co">#&gt; intercept      -4.5813163868 -1.5737283067  5.4453917155  9.934156583</span>
<span class="co">#&gt; dist            2.4086674509  2.7488329575  3.5283145731  4.085425305</span>
<span class="co">#&gt; climb          -0.0093615230 -0.0080203275 -0.0039460163 -0.001159959</span>
<span class="co">#&gt; climb_distance  0.0003751563  0.0005028396  0.0008142368  0.001001573</span>
<span class="co">#&gt; log_sigma2      4.6565098401  4.7253125905  4.9799057245  5.173419493</span>
<span class="co">#&gt;                         75%          95%        97.5%     rhat</span>
<span class="co">#&gt; intercept      14.324050213 20.435459209 23.158888718 1.076772</span>
<span class="co">#&gt; dist            4.639048417  5.520162539  5.830435354 1.078241</span>
<span class="co">#&gt; climb           0.001433777  0.005613192  0.007304871 1.056361</span>
<span class="co">#&gt; climb_distance  0.001203238  0.001502171  0.001603657 1.084557</span>
<span class="co">#&gt; log_sigma2      5.348237447  5.621488633  5.716001010 1.006378</span></pre></body></html></div>
<p>Trace plots provide a visual indication of stationarity. These plots indicate that the MCMC chains are reasonably stationary.</p>
<div class="sourceCode" id="cb8"><html><body><pre class="r"><span class="fu"><a href="../reference/hmclearn-plots.html">mcmc_trace</a></span>(<span class="no">f_hmc</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">2000</span>)</pre></body></html></div>
<p><img src="linear_regression_hmclearn_files/figure-html/unnamed-chunk-8-1.png" width="576"></p>
<p>Histograms of the posterior distribution show that Bayesian parameter estimates align with frequentist estimates for this example.</p>
<div class="sourceCode" id="cb9"><html><body><pre class="r"><span class="fu"><a href="../reference/diagplots.html">diagplots</a></span>(<span class="no">f_hmc</span>, <span class="kw">comparison.theta</span><span class="kw">=</span><span class="no">theta.freq</span>, <span class="kw">burnin</span><span class="kw">=</span><span class="fl">2000</span>)
<span class="co">#&gt; $histogram</span></pre></body></html></div>
<p><img src="linear_regression_hmclearn_files/figure-html/unnamed-chunk-9-1.png" width="576"></p>
</div>
<div id="source" class="section level2">
<h2 class="hasAnchor">
<a href="#source" class="anchor"></a>Source</h2>
<p>A.C. Atkinson (1986) Comment: Aspects of diagnostic regression analysis. <em>Statistical Science</em> 1, 397–402.</p>
</div>
<div id="references" class="section level2">
<h2 class="hasAnchor">
<a href="#references" class="anchor"></a>References</h2>
<p>Venables, W. N. and Ripley, B. D. (2002) <em>Modern Applied Statistics with S</em>. Fourth edition. Springer. ISBN 0-387-95457-0</p>
<p>Gelman, Andrew. (2006) “Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).” <em>Bayesian analysis</em> 1.3: 515-534.</p>
<p>Agresti, A. (2015). <em>Foundations of linear and generalized linear models</em>. John Wiley &amp; Sons. ISBN: 978-1-118-73003-4</p>
<p>Hoffman, M. D., &amp; Gelman, A. (2014). The No-U-Turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. <em>Journal of Machine Learning Research</em>, 15(1), 1593-1623.</p>
</div>
  </div>

  <div class="col-md-3 hidden-xs hidden-sm" id="pkgdown-sidebar">

        <nav id="toc" data-toggle="toc"><h2 data-toc-skip>Contents</h2>
    </nav>
</div>

</div>



      <footer><div class="copyright">
  <p>Developed by Samuel Thomas.</p>
</div>

<div class="pkgdown">
  <p>Site built with <a href="https://pkgdown.r-lib.org/">pkgdown</a> 1.5.1.</p>
</div>

      </footer>
</div>

  


  </body>
</html>
