mathbf x_N\)</span> and response variables <span class="math inline">\(y_1, \dots, y_N\)</span> with likelihood <span class="math display">\[
    p(\mathbf X, \mathbf y | \beta, \beta_0 ) = \prod_{i=1}^N \left[ \frac{1}{1+e^{-\beta_0 + \mathbf x_i \beta}} \right]^{y_i} \left[ 1 - \frac{1}{1+e^{-\beta_0 + \mathbf x_i \beta}} \right]^{1-y_i}
\]</span></p>
<p>First let’s load in the data, we will use the <a href="https://archive.ics.uci.edu/ml/datasets/covertype">cover type dataset</a> commonly used to bench