Stat 550
Daniel J. McDonald
Last modified – 03 April 2024
\[ \DeclareMathOperator*{\argmin}{argmin} \DeclareMathOperator*{\argmax}{argmax} \DeclareMathOperator*{\minimize}{minimize} \DeclareMathOperator*{\maximize}{maximize} \DeclareMathOperator*{\find}{find} \DeclareMathOperator{\st}{subject\,\,to} \newcommand{\E}{E} \newcommand{\Expect}[1]{\E\left[ #1 \right]} \newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]} \newcommand{\Cov}[2]{\mathrm{Cov}\left[#1,\ #2\right]} \newcommand{\given}{\mid} \newcommand{\X}{\mathbf{X}} \newcommand{\x}{\mathbf{x}} \newcommand{\y}{\mathbf{y}} \newcommand{\P}{\mathcal{P}} \newcommand{\R}{\mathbb{R}} \newcommand{\norm}[1]{\left\lVert #1 \right\rVert} \newcommand{\snorm}[1]{\lVert #1 \rVert} \newcommand{\tr}[1]{\mbox{tr}(#1)} \newcommand{\U}{\mathbf{U}} \newcommand{\D}{\mathbf{D}} \newcommand{\V}{\mathbf{V}} \renewcommand{\hat}{\widehat} \]
Carnegie Mellon’s 36-750 Notes
Thank you Alex and Chris for the heavy lifting.
the crash of the Mars Climate Orbiter (1998),
a deadly medical device (1985, 2000),
a massive Northeastern blackout (2003),
the Heartbleed, Goto Fail, Shellshock exploits (2012–2014),
a 15-year-old fMRI analysis software bug that inflated significance levels (2015),
It is easy to write lots of code.
But are we sure it’s doing the right things?
Important
Effective testing tries to help.
(REPL == Read-Eval-Print-Loop, the console, or Jupyter NB)
This tends to result in lots of bugs.
Later on, you forget which values you tried, whether they failed, how you fixed them.
So you make a change and maybe or maybe not try some again.
Write functions.
Lots of them.
👍 Functions are testable
👎 Scripts are not
It’s easy to alter the arguments and see “what happens”
There’s less ability to screw up environments.
I’m going to mainly describe R, but the logic is very similar (if not the syntax) for python, C++, and Julia
fn <- function(formula, data, subset, weights, na.action, method = "qr", model
= TRUE, x = FALSE, y = FALSE, qr = TRUE, singular.ok = TRUE, contrasts =
NULL, offset, ...)fn <- function(e1, e2)fn <- function(.data, ..., .by = NULL, .preserve = FALSE)fn <- function(x, filter, method = c("convolution", "recursive"), sides = 2,
circular = FALSE, init = NULL)fn <- function(n, mean = 0, sd = 1)[1] 0.5855288 0.7094660 -0.1093033[1] 0.5855288 0.7094660 -0.1093033[1] 0.5855288 0.7094660 -0.1093033[1] 0.5855288 0.7094660 -0.1093033hist function
List of 6
$ breaks : num [1:14] -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 ...
$ counts : int [1:13] 4 21 41 83 138 191 191 182 74 43 ...
$ density : num [1:13] 0.008 0.042 0.082 0.166 0.276 0.382 0.382 0.364 0.148 0.086 ...
$ mids : num [1:13] -2.75 -2.25 -1.75 -1.25 -0.75 -0.25 0.25 0.75 1.25 1.75 ...
$ xname : chr "rnorm(1000)"
$ equidist: logi TRUE
- attr(*, "class")= chr "histogram"[1] "histogram"[1] 3[1] 2 3 4 5Error in x + 1: non-numeric argument to binary operatorincrementer <- function(x, inc_by = 1) {
if (!is.numeric(x)) {
stop("`x` must be numeric")
}
x + 1
}
incrementer("a")Error in incrementer("a"): `x` must be numeric[1] 3incrementer <- function(x, inc_by = 1) {
if (!is.numeric(x)) {
stop("`x` must be numeric")
}
x + inc_by
}
incrementer(2, -3)[1] -1library(testthat)
incrementer <- function(x, inc_by = 1) {
if (!is.numeric(x)) stop("`x` must be numeric")
if (!is.numeric(inc_by)) stop("`inc_by` must be numeric")
x + inc_by
}
test_that("incrementer validates arguments", {
expect_error(incrementer("a"))
expect_equal(incrementer(1:3), 2:4)
expect_equal(incrementer(2, -3), -1)
expect_error(incrementer(1, "b"))
expect_identical(incrementer(1:3), 2:4)
})── Failure: incrementer validates arguments ────────────────────────────────────
incrementer(1:3) not identical to 2:4.
Objects equal but not identicalError:
! Test failedA unit is a small bit of code (function, class, module, group of classes)
A test calls the unit with a set of inputs, and checks if we get the expected output.
gcd <- function(x, na.rm = FALSE) {
if (na.rm) x <- x[!is.na(x)]
if (anyNA(x)) return(NA)
stopifnot(is.numeric(x))
if (!rlang::is_integerish(x)) cli_abort("`x` must contain only integers.")
if (length(x) == 1L) return(as.integer(x))
x <- x[x != 0]
compute_gcd(x) # dispatch to a C++ function
}
test_that("gcd works", {
# corner cases
expect_identical(gcd(c(1, NA)), NA)
expect_identical(gcd(c(1, NA), TRUE), 1L)
expect_identical(gcd(c(1, 2, 4)), 1L)
# error
expect_error(gcd(1.3))
# function
expect_identical(gcd(c(2, 4, 6)), 2L)
expect_identical(gcd(c(2, 3, 7)), 1L)
})Unit testing consists of writing tests that are
Among others:
Each test focuses on a single component.
Named so that you know what it’s doing.
Collection of related tests in a common context.
Prepares the environment, cleans up after
(loads some data, connects to database, necessary library,…)
Test suites are run and the results reported, particularly failures, in a easy to parse and economical style.
For example, Python’s {unittest} can report like this
$ python test/trees_test.py -v
test_crime_counts (__main__.DataTreeTest)
Ensure Ks are consistent with num_points. ... ok
test_indices_sorted (__main__.DataTreeTest)
Ensure all node indices are sorted in increasing order. ... ok
test_no_bbox_overlap (__main__.DataTreeTest)
Check that child bounding boxes do not overlap. ... ok
test_node_counts (__main__.DataTreeTest)
Ensure that each node's point count is accurate. ... ok
test_oversized_leaf (__main__.DataTreeTest)
Don't recurse infinitely on duplicate points. ... ok
test_split_parity (__main__.DataTreeTest)
Check that each tree level has the right split axis. ... ok
test_trange_contained (__main__.DataTreeTest)
Check that child tranges are contained in parent tranges. ... ok
test_no_bbox_overlap (__main__.QueryTreeTest)
Check that child bounding boxes do not overlap. ... ok
test_node_counts (__main__.QueryTreeTest)
Ensure that each node's point count is accurate. ... ok
test_oversized_leaf (__main__.QueryTreeTest)
Don't recurse infinitely on duplicate points. ... ok
test_split_parity (__main__.QueryTreeTest)
Check that each tree level has the right split axis. ... ok
test_trange_contained (__main__.QueryTreeTest)
Check that child tranges are contained in parent tranges. ... ok
---------------------------------------------------------
Ran 12 tests in 23.932sR exampleCore Principle:
Tests should be passed by a correct function, but not by an incorrect function.
The tests must apply pressure to know if things break.
Make sure that incorrect functions won’t pass (or at least, won’t pass them all).
Tip
Cover all branches.
Make sure there aren’t branches you don’t expect.
Assertions are things that must be true. Failure means “Quit”.
def fit(data, ...):
for it in range(max_iterations):
# iterative fitting code here
...
# Plausibility check
assert np.all(alpha >= 0), "negative alpha"
assert np.all(theta >= 0), "negative theta"
assert omega > 0, "Nonpositive omega"
assert eta2 > 0, "Nonpositive eta2"
assert sigma2 > 0, "Nonpositive sigma2"
...The parameters have to be positive. Negative is impossible. No way to recover.
Errors are for unexpected conditions that could be handled by the calling code.
Code can also do this. It can try the function and catch errors to recover automatically.
For example:
Load some data from the internet. If the file doesn’t exist, create some.
Run some iterative algorithm. If we haven’t converged, restart from another place.
Code can fix errors without user input. It can’t fix assertions.
Test Driven Development (TDD) uses a short development cycle for each new feature or component:
Write tests that specify the component’s desired behavior.
The tests will initially fail because the component does not yet exist.
Create the minimal implementation that passes the test.
Refactor the code to meet design standards, running the tests with each change to ensure correctness.
Writing the tests may help you realize
The tests define a specific plan for what the function must do.
You will catch bugs at the beginning instead of at the end (or never).
Testing is part of design, instead of a lame afterthought you dread doing.
test_1 doesn’t help you at all, but test_tree_insert makes it easy for you to remember what the test is for.
Do this
foo <- function(x) {
if (x < 0) stop(x, " is not positive")
}
foo <- function(x) {
if (x < 0) message(x, " is not positive")
# not useful unless we fix it too...
}
foo <- function(x) {
if (x < 0) warning(x, " is not positive")
# not useful unless we fix it too...
}
foo <- function(x) {
if (length(x) == 0)
rlang::abort("no data", class="no_input_data")
}These allow error handling.
See here for more details.
Seems like overkill,
but when you run a big simulation that takes 2 weeks,
you don’t want it to die after 10 days.
You want it to recover.
tib <- tibble(
x1 = rnorm(100),
x2 = rnorm(100),
y = x1 + 2 * x2 + rnorm(100)
)
mdl <- lm(y ~ ., data = tib)
class(tib)[1] "tbl_df" "tbl" "data.frame"[1] "lm"The class allows for the use of “methods”
R “knows what to do” when you print() an object of class "lm".
print() is called a “generic” function.
You can create “methods” that get dispatched.
For any generic, R looks for a method for the class.
If available, it calls that function.
There are lots of generic functions in R
Common ones are print(), summary(), and plot().
Also, lots of important statistical modelling concepts: residuals() coef()
(In python, these work the opposite way: obj.residuals. The dot after the object accesses methods defined for that type of object. But the dispatch behaviour is less robust.)
The convention is that the specialized function is named method.class(), e.g., summary.lm().
If no specialized function is defined, R will try to use method.default().
For this reason, R programmers try to avoid . in names of functions or objects.
Call:
lm(formula = y ~ ., data = tib)
Coefficients:
(Intercept) x1 x2
0.1216 1.0803 2.1038 [1] "This is an linear model."The advantage is that you don’t have to learn a totally new syntax to grab residuals or plot things
You just use residuals(mdl) whether mdl has class lm could have been done two centuries ago, or a Batrachian Emphasis Machine which won’t be invented for another five years.
The one draw-back is the help pages for the generic methods tend to be pretty vague
Compare ?summary with ?summary.lm.
These are often tricky, but are very common.
Most programming languages have this concept in one way or another.
In R code run in the Console produces objects in the “Global environment”
You can see what you create in the “Environment” tab.
But there’s lots of other stuff.
Many packages are automatically loaded at startup, so you have access to the functions and data inside those package Environments.
For example mean(), lm(), plot(), iris (technically iris is lazy-loaded, meaning it’s not in memory until you call it, but it is available)
Other packages require you to load them with library(pkg) before their functions are available.
But, you can call those functions by prefixing the package name ggplot2::ggplot().
You can also access functions that the package developer didn’t “export” for use with ::: like dplyr:::as_across_fn_call()
As one might expect, functions create an environment inside the function.
Non-trivial cases are data-masking environments.
Error in eval(expr, envir, enclos): object 'x2' not foundlm() looks “inside” the tib to find y and x2lm() environmentWhen Knit, .Rmd files run in their OWN environment.
They are run from top to bottom, with code chunks depending on previous
This makes them reproducible.
Jupyter notebooks don’t do this. 😱
Objects in your local environment are not available in the .Rmd
Objects in the .Rmd are not available locally.
Tip
The most frequent error I see is:
The reason is almost always that the chunks refer to objects in the Global Environment that don’t exist in the .Rmd
library() calls were made globally but not in the .Rmd
.Rmd in your file system
?function to see the help
.Rmd won’t Knit
rstudioapi::restartSession(), then run the Chunks 1-by-1Adding browser()
browser() to the code somewherebrowser() and you’ll have access to the local environment to play aroundQuestion I get uncountably often that I hate:
“I ran the code like you had on Slide 39, but it didn’t work.”
Don’t just share a screenshot!
Unless you get lucky, I won’t be able to figure it out from that.
And I can’t copy-paste into Google.
What you need is a Reproducible Example or reprex. This is a small chunk of code that
The best way to do this is with the {reprex} package.
Open a new .R script.
Paste your buggy code in the file (no need to save)
Edit your code to make sure it’s “enough to produce the error” and nothing more. (By rerunning the code a few times.)
Copy your code.
Call reprex::reprex() from the console. This will run your code in a new environment and show the result in the Viewer tab. Does it create the error you expect?
If it creates other errors, that may be the problem. You may fix the bug on your own!
If it doesn’t have errors, then your global environment is Farblunget.
The Output is now on your clipboard. Share that.
Note
Because Reprex runs in it’s own environment, it doesn’t have access to any of the libraries you loaded or the stuff in your global environment. You’ll have to load these things in the script. That’s the point
Suppose we want to find \(\max_x f(x)\).
We repeat the update \(x \leftarrow x + \gamma f'(x)\) until convergence, for some \(\gamma > 0\).
Goal: find the MLE for \(\lambda\) using gradient ascent
\[ \begin{aligned} L(\lambda; y_1,\ldots,y_n) &= \prod_{i=1}^n \frac{\lambda^{y_i} \exp(-\lambda)}{y_i!}I(y_i \in 0,1,\ldots)\\ x &\leftarrow x + \gamma f'(x) \end{aligned} \]
UBC Stat 550 - 2024