R the Software Engineering Way: Chapter 2

We're up to Chapter two already, and I can't think of a better way to start the new year! As always, comments, feedback and requests for help are more than welcome. In the interests of paying my rent and my bills, I'm also on the lookout for contracting and consulting roles: given this, ones involving R, DevOps, containerisation or sysadmin things with Linux might be good fits. The usual donation channels are also available, as always, should you wish to help.

Chapter Zero

Chapter One

In the last two chapters of this book we've gone through the basics of setting up a development container, building a package and writing some functions that perform basic mathematical tasks. What we've written, however, is all pretty trivial. Now, while teaching how to properly work with a truly large codebase is something that would take a book to teach in itself, we can certainly build something a bit more complex than the toy problems we've worked with so far. The question is: what?

The commit for this chapter can, as previously, be found at https://gitlab.com/irishenceaway/r.the.software.way: the commit for this chapter is tagged chapter2.

The problem

What problem to choose for teaching purposes is not a question easily answered. After all, there are multiple desiderata, and some of them are in conflict with each other. For starters, we want something that actually requires you to write code rather than simply importing a package in its entirety: you need to actually use the control flow elements that we discussed last chapter. We also want something that has actual bearing on statistical methods that you might use day-to-day: while I could make you implement a Finite Difference scheme or something, the average statistician probably isn't very interested in Partial Differential Equations. Finally, whatever is implemented can't be too complex, because I'm not trying to teach statistics, really, and even in the context of software engineering we don't want to be writing truly complex logic in a learning situation (handling that is something you can only learn through bitter experience).

Fortunately there's a statistical method that ticks all of these boxes: Approximate Bayesian Computation. While there are packages available for it, they're never quite as easy to use as the corresponding packages for frequentist statistics, and even the best packages available require more manual work than a frequentist method will, because Bayesian methods require quite a lot of individual judgement and decision-making in choice of priors, cutoffs and suchlike. Even when a method is available, then, it's often easier to write simple ABC schemes yourself in a lot of cases, or to write parts of it yourself while using packages for other parts. Consequently, what you'll write here isn't a toy problem: it's a simple implementation, but it's very much something that you could, if faced with a situation where you need to, build on to use in a real application. In this chapter, then, we'll be implementing a basic rejection sampling scheme for Approximate Bayesian Computation.

While this may sound a little intimidating, we have built up all the skills that we need to write an implementation of this over the last few chapters: we have the package management skills, the understanding of functions and the control flow tools that we need to do this. And on the subject of package management...

A better package manager: using rv

At this point, we encounter an issue, though it might not be immediately obvious: while the last two chapters only really used the devtools package, which is unusual in a bunch of ways (it really is more of a cli utility masquerading as a package than a package in itself), in this chapter we're going to have to import quite a few actual packages that we will use for actual R tasks. These packages will have versions that need to be tracked and will have a whole bunch of dependencies, which means that our current model for handling packages is no longer really going to cut it. Also, for better or worse (I personally don't like it), we do a lot of stuff in R from within the interactive shell. This means that even when we're doing almost all of our work in the development container we still find ourselves maintaining a local environment, which, if not isolated, can lead to some weird behaviour.

What I really wanted was a tool like uv (a python package manager) or npm (a project-based package manager for JavaScript). These package manager are run from the command line (and thus mean that we can manage dependencies without having to start up an interactive session, which works more nicely with cli-based workflows than what R usually does), are project-based, so we get project-level package isolation, and unlike basically any other R package manager I've seen, are declarative (so you always wind up with a working set of project dependencies if it's possible). Unfortunately, I didn't think any such package manager existed. But then lo, the heavens opened and revealed to me the existence of the perfect package manager for my purposes: rv ¹.

rv is a package manager that's quite different from basically every other R package manager available. For a start, while other package managers are called from within an interactive R shell (we saw this with both pak and the default install.packages functionality), rv runs from the command line in the same way that tools like npm and uv do. While a statistician might find this workflow a bit annoying initially, it is a much better and more usable workflow from the perspective of a software engineer.

Secondly, and more importantly, rv is declarative rather than being imperative. Most R package managers install packages one-by-one, with dependencies being installed each time a package asks for them to be installed. The package manager then retroactively creates a snapshot showing what the state ended up being after everything's done. While this works, more or less, it can lead to some strange outcomes, such as one package overwriting another when it installs a dependency, causing the final state to end up broken. It can also make replicability difficult, as different systems can plausibly install things in different orders, leading to different end states. Taken together, this means that even assiduous package management will not guarantee replicability.

The declarative approach, by contrast, starts with describing a desired end state: "I want this set of packages installed, in these versions, handling dependencies in this way." rv then analyses the state that you've written and figures out a plan for getting you there. If getting to the end state is impossible, rv will tell you rather than, as is often the case, breaking your install. This leads, in the end, to much better replicability and fewer headaches all around.

There are a bunch of other benefits: when you open an R session in an active rv project, your packages will often be loaded automatically (to an extent: this doesn't always work as you might expect). rv is also gratifyingly fast: while CRAN downloads can still take some time, the rv cache is pretty efficient and the tool itself is written in Rust, making it pretty easy to use. Finally, rv is by far the best option available if, like us, we do most of our work in a development container but also want to maintain a local environment that's similar to the containerised one. All in all, it's just a much more powerful tool than the others that are available.

To install rv locally, you'll simply run:

curl -sSL https://raw.githubusercontent.com/A2-ai/rv/refs/heads/main/scripts/install.sh | bash

In your command line: this uses cURL to pull the install script from a Github repository and pipes the script to bash (remember pipes from Chapter 0? This is one of the things they're good for). This isn't a very secure installation method in a bunch of ways, as you have to trust the URL not to execute arbitrary stuff on your system, but needs must, and this particular install seems safe at the time of writing. After the install you'll need to restart the shell (more on this later), and then you can begin to use rv.

rv, like uv or renv, is a project-based package manager and requires you to create a project before you can start using it. To do so, you'll want to run rv init . in your project directory to initialise it as an rv project (this will quickly get weird if you're already using a package manager, so maybe don't: if you've used renv, though, you can run rv migrate renv to migrate the project effectively). Doing this will create two files called rproject.toml and rv.lock respectively. Now, in many cases rv won't configure what repositories to draw from initially, so you'll be unable to install packages until you do that. To do that, you might start by opening up the rproject.toml file. The contents will look like this:

[project]
name = "r.the.software.way"
r_version = "4.5"

# A list of repositories to fetch packages from. Order matters: we will try to get a package from each repository in order.
# The alias is only used in this file if you want to specifically require a dependency to come from a certain repository.
# Example: { alias = "PPM", url = "https://packagemanager.posit.co/cran/latest" },
repositories = [
]

# A list of packages to install and any additional configuration
# Examples:
    # "dplyr",
    # {name = "dplyr", repository = "CRAN"},
    # {name = "dplyr", git = "https://github.com/tidyverse/dplyr.git", tag = "v1.1.4"},
dependencies = [
]

You will note the repositories entry, which is currently an empty list. This means that we currently have no repositories to install packages from. To change this, you can change the repositories entry to look like this:

repositories = [
    { alias = "PPM", url = "https://packagemanager.posit.co/cran/latest" },
]

This adds the Posit package manager repository as a package source to rv: this is probably the best option for packages available at the moment, and hopefully it will remain that way for the foreseeable future if Hadley Wickham keeps doing good work. You can now run rv add <package-name> (from the command line, not from within an R session) to install a package. You can then proceed to start an interactive R session, and, mirabile dictu, you will find the package available to you (you may or may not have to install it). This allows you to run a more-or less fully reproducible environment for your R scripts locally, which, as we've established, is quite useful for a lot of short and otherwise ad-hoc tasks or for using tooling which in most environments would live in the command line but that R does interactively.

However, one way or another we still do most of our work in the development container, so getting rv set up and reproducing the project in the container is our next goal. And this is a bit more of a challenge (though I will stress that it's the good kind of challenge that will teach you about the command line, Linux and container technology, not the kind that makes you want to tear your hair out).

Getting your container image working with rv

As I have alluded to immediately earlier, getting rv running in a container can take a bit of doing. We can start by looking at the development image's Dockerfile post-rv:

FROM rocker/r-ver:4.5.2 AS base

RUN apt-get update

RUN apt-get install -y curl libcurl4-openssl-dev libicu-dev libtiff-dev make zlib1g-dev libssl-dev libx11-dev libxml2-dev libfontconfig1-dev libharfbuzz-dev libfribidi-dev libfreetype6-dev libpng-dev libtiff5-dev libjpeg-dev libwebp-dev pandoc

RUN curl -sSL https://raw.githubusercontent.com/A2-ai/rv/refs/heads/main/scripts/install.sh | bash

RUN mv ~/.local/bin/rv /usr/local/bin/rv

WORKDIR /r.the.software.way

COPY rproject.toml rproject.toml

COPY rv.lock rv.lock

RUN rv sync

RUN rv activate

A few things have become simpler: we're no longer calling any R at all during the build process and we can dispense with the packages.R file entirely. That said, things have otherwise gotten quite a bit more complex. First off, the systems packages: the tidyverse and devtools rely on quite a lot of systems packages that we have to install via the systems package manager. We did this for libcurl in the last Dockerfile, but here we have to install quite a few more, because while with the tools we were using previously the systems libraries weren't necessary for things to work (the packages install, but fail silently), rv, being declarative in nature, will actively fail to install the package unless all systems dependencies are present. This may initially seem like a pain in the behind, but it saves a lot of time and pain when debugging down the line.

Next, we install the rv executable in the usual way. There's a catch here, though: remember how you had to restart the terminal to get the tool to work locally? That's because in order to find command line executables, the shell that you're working with looks at an environment variable called PATH. You can have a look at your path variable by executing

echo $PATH

in a terminal. Now, here's the thing: rv installs itself in a local user directory and adds itself to the path using the .bashrc file, which is a file that contains commands that run whenever you open a new shell (it's essentially your profile configuration). On a development machine, that usually works fine, but in a container image, a lot of the infrastructure that rv takes for granted isn't present and the PATH variable often doesn't act as you might expect. We could try and alter PATH to make it work, but the easier option is simply to copy the executable (as rv is written in Rust, this is a single binary file) to where all of the other executables live using

RUN mv ~/.local/bin/rv /usr/local/bin/rv

This copies the executable from the place where rv installed it to /usr/local/bin/rv, which is one of the locations where executables can live on Unix-like systems.

Finally, we can copy across rproject.toml and rv.lock and then run rv sync to install exactly the same versions of exactly the same packages that you have installed locally into your container. We then run rv activate to activate the project, and we're good to go!

Before we use our new image, we'll need to update our compose.yaml a little. Previously we just bound our entire project directory to the container, and while this works, things can get a little weird if we bind the local versions of rproject.toml and rv.lock over the container versions. We thus need to update our compose.yaml file to look more like this:

services:
  r-development-container:
    image: iris-meredith/r-the-software-way
    stdin_open: true
    tty: true
    volumes:
    - type: bind
      source: R
      target: /r.the.software.way/R
    - type: bind
      source: tests
      target: /r.the.software.way/tests
    - type: bind
      source: data
      target: /r.the.software.way/data
    - type: bind
      source: scripts
      target: /r.the.software.way/scripts
    - type: bind
      source: plots
      target: /r.the.software.way/plots
    - type: bind
      source: DESCRIPTION
      target: /r.the.software.way/DESCRIPTION
    - type: bind
      source: NAMESPACE
      target: /r.the.software.way/NAMESPACE

Here, rather than binding everything, Docker Compose synchronises only the directories and files that we've specified, ignoring the rv library and the .vscode directory. This avoids binding over the library directories or any files that rv relies on for managing its project, while still binding the source code and tests so that we use the development container as usual. We can then start the container in the usual way by running

sudo docker compose run --rm r-development-container

and find ourselves in an interactive R shell with exactly the same environment that we have locally, all neatly packaged and reproducible!

Exercise: get set up with rv

In your project directory:

• Install rv and make sure that it's on your path and that you can run it.
• Run rv init in your project directory to initialise an rv project.
• Install the tidyverse and devtools packages.
• Run the existing tests locally to confirm that everything works.
• Modify your Dockerfile and compose.yaml files to install the packages from rproject.toml into your image and use them properly.
• Run the tests again in the development container to confirm that everything works as intended.
• Commit your changes to the repository and push them.

Handling data, scripts and plots

As we've noted earlier, almost all of the R code for a package lives in the /R directory as function declarations. This is mostly fine: we can write the function definitions, and even if we haven't built the package we can use devtools to import our functions and call them as normal from an interactive session: not a hugely difficult task, all told.

However, we'll probably, in the course of an analysis, want to write some scripts that we don't want to distribute with the package (should we choose to put it on CRAN or something). This means that we need somewhere to put them. A similar issue arises with plots, which you'll probably create a few of in the exploratory analysis step. Quite the opposite issue arises with data: you'll often want to distribute data along with your package.

Fortunately, the solution to this is quite simple: any data you want to distribute with the package lives in /data, and for files that we don't want to distribute with the package we can create /scripts and /plot directories which the R packaging system will ignore. We can bind those to the container as usual, which will have the additional benefit of making any generated plots or files available to us outside of the container: useful if you want to, for example, insert the plots into a report or a slideshow.

Git interlude: branching, merging and .gitignore

By now you'll most likely have gotten quite used to the git add, git commit, git push workflow that we've used in the last few chapters. On the opposite side of things, you've likely also become quite familiar with the git clone command. This is, for the most part, all that you really need to know, but with the work we're doing becoming more complex, some more complex uses of git might be worth discussing. This section thus discusses two things of some importance: the branch/merge workflow and the .gitignore and .gitkeep files.

.gitignore and .gitkeep

As we've discussed, this analysis is likely to require you to make some plots. Plots are important for most statistical work, in fact, and you're likely to want to make and save a lot of them. However, plots don't version control well: they're big binary files that can take a while to upload and download, you can't really track changes in them in any meaningful way and they're just generally a pain. Besides, if you have the code to generate them stored in /scripts or somewhere in your package functions, they're trivial to recreate. We thus want to tell git not to track any image files in question, and we can do this with a .gitignore file. This file is called .gitignore and lives in the root directory, and it looks something like this:

*.png
*.jpg
*.jpeg
*.webp
*.pdf

This example .gitignore file will tell git not to track any images that you create, as well as any PDFs that you might generate. While in this case we're simply trying to make things a bit faster and better-behaved, .gitignore is often also used to prevent things like .env files, that often contain secrets like API keys that you want to keep out of public repositories, from being tracked by git.

Another important git file is the .gitkeep file. When we created the /scripts and /plots directories, they start off empty until you do something with them. In the normal course of events, this means that they won't be tracked by git. If you want git to track those directories, then, you can create an empty .gitkeep file in the directory to force git to track it.

Branching and merging

Thus far, we've only really worked on one branch, which we've called main. All of our commits and pushes have been on this branch. For small projects or ones maintained by a single person, this generally works. However, if you're working on multiple things at once (common enough when you're working on a new feature, have to stop to fix a bug in your current release and then want to go back to working on the feature) or if you're working on a codebase with multiple people, it's worth understanding the basics of branching and merging in git.

Those of you who've been pulling the accompanying repository to follow along with this book will already have created new branches once or twice: the git checkout <tag> -b <branch-name> command that you use to check out a version of the repository corresponding to a chapter tag creates a new branch with a name that you've given it and switches to it. A branch in git is simply a pointer to a given snapshot of the repository that we've given a name to and can now reference. git checkout, however, is a bit of a dated command (I've only used it because I don't want to think too hard about how git switch works with tags). Here, then, is how we might execute a branching workflow in modern git.

We start off on the main branch. Now, let us assume that we're working through the exercises in this chapter, and we don't want to do this on the main branch (maybe it's important to us that we're able to pull and use that code as is without half-finished chapter 2 exercises). We thus run git switch -c chapter2, which creates a new branch called chapter2. We do some work on the exercises, and then find a bug in the chapter 1 work which has to be fixed immediately. We don't want to lose our in-progress work on chapter 2, so we'll either want to commit our current changes with git commit, or if we aren't ready to commit yet, we can run git stash to save uncommitted changes and rollback to the last commit. Either way, we can then run git switch main to switch back to the main branch, which contains your code as it stood at the end of chapter 1. We can then create another branch by running git switch bugfix or something more descriptive (definitely something more descriptive, really). We write whatever code we need to fix the bug, then run git commit. Now, this commit currently exists only on the bugfix branch, so we'll need to merge it back into main. To do this, we run git switch main followed by git merge bugfix: this will merge the bug fix from the bugfix branch into main. We can then push main to our remote repository and resume our work on chapter 2 by running git switch chapter2, followed by git stash pop to restore the uncommitted work.

Once you've finished work on chapter2, you'll want to merge it back into main, and at this point, a potential issue can arise. If the bugfix and the chapter 2 work were in different files, or in the same file but relatively simple, there won't be issues: however, if both branches have modified the same code in different ways, you might find yourself with a merge conflict on your hands: the same lines have been changed in two different ways. In this case, you'll usually be prompted to resolve the merge conflict manually, by figuring out what code captures the behaviour you want and setting it to that state. Merge conflicts aren't usually too thorny to sort out: they can get ugly in big codebases and people have horror stories about them, but in 95% of cases they're easily resolved. If you want to check for possible merge conflicts before trying to merge, the git diff command can be quite helpful.

Finally, after a period of working like this, you'll often wind up with a lot of branches that you no longer use. While retaining the chapter branches might be helpful (though of course we have tags for that), we probably don't want to retain every out-of-date branch for every bugfix we've ever performed. To list existing branches, then, we can run git branch, and to delete an unused branch we can run git branch -d <branch-name>. For example, deleting the bugfix branch would use git branch -d bugfix.

As for remote repositories, you can push other branches to your remote using the usual git push origin <branch-name> command. There are a number of reasons to do this: you might, for example, wish to maintain a dev branch separate to your main branch where you keep code that's been less well tested and reviewed than your stuff on main. For now, though, keeping all of your branches except for main local will probably work well enough.

While this workflow can be a bit of a pain to work with at first (not much, but a little), in the long run this will make managing workflows a lot easier and prepare you effectively for working on a codebase with other people. It's thus highly worth mastering. And now, with all of that in place, we can finally get around to actually writing the rejection sampler!

Writing the rejection sampler

The statistical underpinnings

The basic principle behind ABC is simple. Suppose, to begin with, that we have a vector of data that we think we can fit some statistical model to. This model will have some set of parameters that are unknown to us and that we want to predict: for example, if we think the data is generated by a Gaussian process, our model will have two parameters: . In most statistical approaches, we would use a method here that relies on the likelihood function of the data , or at least (as in the case of Markov Chain Monte Carlo methods), a function proportional to the likelihood. However, a closed-form function for this isn't always available, or if available, may be intractable in several ways. In this case, with your average statistical methods, you're kinda fucked.

Approximate Bayesian Computation takes as its starting point the fact that even when you don't have a closed-form function for the likelihood and are unable to get one, you can often generate simulated data from a hypothetical process with the given parameter. Given this, it's possible to perform rejection sampling: we generate random values for the parameters from our prior distributions and simulate a dataset from said parameters. We then calculate a distance . If , where is a small parameter that you set when building the ABC sampler, we add to our dataset of estimates of ; otherwise, we reject it. We repeat this process until we have a specified number of estimated values: the resulting dataset is a sample from the posterior distribution, and in practice we'd usually take the mean or median of that sample to be our parameter estimate.

This approach is robust and powerful, but there are some challenges that will need to be considered in the code. First off, this is a Bayesian method, which means that when doing parameter estimation, we need to choose priors (I apologise sincerely if you're coming at this from a frequentist background and I've just spooked you: I assure you it's not too hard). The choice of priors can both bias your results (though in the limit as the number of samples tends to infinity it should converge to the correct value unless you've chosen something really fucked up), and, more worryingly, make getting your sample very slow: after all, if you find yourself sampling from a prior that is completely different from what the parameters actually are (or that possibly doesn't even include the actual parameter value at all), it will take you a very long time to get the number of samples that you want. Secondly, the choice of distance metric or divergence is left entirely up to the analyst, and depending on your exact problem, there are a whole lot of different things you could try. One popular approach is to calculate summary statistics and find the distance between them, but that merely defers the problem. Another approach (and my preferred one) is to use a metric that avoid the issue: the Wasserstein distance metric is pretty good for this most of the time, but can get expensive for high dimensions. Of course, you also have to choose the cutoff parameter, and that's in itself a delicate balancing act between choosing one small enough to give you good results and choosing one large enough that you get the samples you want in a reasonable time.

While many of these issues are primarily statistical, these also need to be addressed in the design of the sampler. It is important, for example, to be able to decouple your priors and your simulation code from the core loop, as you're likely to change them multiple times during the analysis and having to edit the core loop ad-hoc each time is going to cause headaches. You'll also want an easy way to adjust your parameters and the number of parameter samples you want. Finally, there are very real questions of performance here: if you write something inefficiently you can find yourself with a process eating multiple CPU cores while doing essentially nothing.

Writing the code

With this in mind, we'll want to break down the ABC rejection sampling scheme into the following steps:

• Select a model for your data: you'll usually do this after performing some basic exploratory analysis like plotting a histogram, running heteroscedasticity tests and computing some basic summary statistics. This usually comes down to the choice of the statistician.
• With the model in hand, select priors for your parameters. These should reflect information that you have about your data: for example, variances can't be negative, so when choosing a prior distribution for your variance parameter, you'll want to draw from a distribution that's strictly non-negative. There is, however, no theoretical upper bound to a variance, so something like a fairly flat exponential distribution might be appropriate if you're looking for a maximally uninformative prior, or you can go with a tighter one if you have more of an idea as to what a plausible value for the variance might be.
• Write a function that takes as arguments your parameters and generates a simulated dataset of a given length.
• Choose or write your distance metric. I've already talked a bit about options for that: your exact choice will depend on your problem.
• Write a function that takes as arguments your data, your prior function, your simulator, your distance metric (which should also be a function), your cutoff and the number of parameter estimates you wish to generate. This function should have a while loop as its core that contains the accept/reject logic described above and returns an object (this being R, probably a tibble) containing the accepted parameter estimates.
• Finally, you'll have written all this code so far as functions in your package. This makes it a bit irritating to run, so you'll want to create a short script in the scripts file that imports your functions and runs the rejection sampler with your choice of parameters. An alternative approach is to write a main function in your package, but I personally don't like this solution.

This breakdown has divided the overall problem into neat, isolated functions that we can write by themselves with limited reference to other parts of the program. Obviously the function interfaces have to match up, but the internals can be written safely without having to think about anything beyond the interface, making all of this testable, modular and easier to modify than it otherwise would be.

Testing the code

Seeds and snapshot tests

When writing tests for the code you've written, you'll quickly run into an obvious difficulty: unlike the code we wrote in Chapter 1, much of this code is non-deterministic in output.

Fortunately testthat provides one part of the solution to the puzzle in the form of snapshot tests, and withr provides the other in the form of local_seed.

As you might know, the random samples generated by R's distribution functions aren't truly random: rather, they're what we call pseudorandom. This means that they have almost all of the statistical properties of random numbers (though if you feel like reading up about the differences, One-Time Pad encryption is both fascinating to learn about and is an important case where pseudorandomness isn't good enough), but they're generated through a deterministic software process. In particular, the pseudorandom number generators rely on an initial value called a seed: if a pseudorandom number generator is called twice with the same seed set, it will generate the same sequence of outputs each time. R usually uses the system clock as its seed: this ensures that every run of a random number generator produces a different output. However, for testing purposes, we can set a seed so that we get reproducible results each time we run them.

The simplest way to do this is to set the seed globally before running tests with set.seed. However, this leaves something to be desired: if you forget to unset the seed afterwards, anything else that you do in the global session might behave strangely. What you can do instead is use the withr library to set a seed for a particular codeblock:

with_seed(0, rnorm(1,0,1))

This, when executed, will return exactly the same value every time, without changing the global seed: if you execute the normal distribution sampler outside of this block, it will generate values nondeterministically as usual. This, then, neatly solves the random number problem in testing.

The second issue, though, is that manually constructing test cases for this manner of problem can be a real headache: it is, after all, simply not practical to manually generate a simulation output, then manually calculate whatever distance metric you chose, then manually... you get the point, anyway. In order to avoid this, we can use testthat's snapshot test feature. For example, if we wanted to test the fizzbuzz_list function that we wrote in Chapter 1 using a snapshot test, we can add an extra line to the test-fizzbuzz-list.R file as follows:

expect_snapshot(fizzbuzz_list(100))

The first time this test is executed, testthat will create a markdown file containing the output of this file as reference output in tests/testthat/_snaps. During future runs, the tests will compare the output against the expected output, passing the test if it's identical and failing it if it isn't. Notably, this allows for much larger volumes of output to be tested: the equality tests we wrote last chapter only tested a couple of individual elements, while the snapshot tests can test that the entire output for a full hundred numbers is as it should be. The catch with snapshot tests, though, is that they rely on the code you've written being right the first time you've written it: it's thus worth, after running the tests for the first time, to review the snapshot file in detail and confirm that the output is correct for the given arguments.

Using these two tools together, we can now write a test for our rejection sampling code. Suppose that we have a function called sample_priors_mixture (which you will shortly write for yourself) that generates sample parameters from the joint prior distribution. A test for that function might look like this:

test_that("Mixture-of-Gaussian priors are sampled correctly", {
  library(withr)
  with_seed(0, expect_snapshot(sample_priors_mixture()))
})

In this case, we've wrapped the expect_snapshot statement inside a with_seed statement that guarantees that the snapshot test runs with a seed of zero each time, then executes the snapshot test in the usual way. The contents of the snapshot file for this test looks like the following:

# Mixture-of-Gaussian priors are sampled correctly

    Code
      sample_priors_mixture()
    Output
              x1         x2     sigma1     sigma2 proportion
       7.9339440 -4.6898267  2.3632856  0.2914135  0.9082078

This is correct for the function we've written so we can accept it and leave it as is. Running the tests again, we can see that this replicates correctly.

As you grow more experienced with testing, you'll develop a sense for when you need to write tests and when you don't: in general if the code is very well-understood and simple, you shouldn't spend resources on testing it. Until you develop that judgement, however, I highly recommend testing everything just to help you get into the habit of doing it.

Exercise: write a basic rejection sampler

In this exercise, you'll be working with two sets of data that can be found in the /data folder of the repository: abc_gaussian_data.csv and abc_mixture_data.csv. The first data is a sample from a Gaussian distribution, and the second is a sample from a mixture of two Gaussians. Assume, for the sake of this exercise, that finding closed-form likelihoods for these problems is impossible.

Your goal is to write a set of functions that will allow you to estimate the parameters of the distributions that these datasets were drawn from. To do this, we'll follow the steps discussed above:

• Create a new development branch of your repository and switch to it. Stash any uncommitted changes first.
• In a script or an interactive session, perform basic exploratory analysis on the abc_gaussian_data.csv dataset. Generate some plots and store them in the correct directory.
• On the basis of these plots, select a joint prior distribution for parameters and and write a function to generate a sample from it. Write a unit test to test this and run the test to generate a snapshot. Run the tests again and confirm that the output is deterministic.
• Write a function to generate simulated data from a Gaussian distribution with the given prior parameters. Make sure that the interface matches the previous function.
• Write a function to calculate a suitable distance metric between the real sample and the simulated sample, or find a package that makes one available. If you use a package, install it using rv and rebuild your container image. If you write it yourself, write a unit test for it.
• Write a function that takes the three previous functions as arguments, along with whatever else you need, and implements the rejection sampling scheme described above. It should return a tibble containing your parameter estimates.
• Find a set of hyperparameters that gives you decent convergence and decently accurate estimates of your parameters. Once you've done this, write a snapshot test for the function and verify it as before.
• Plot your parameter estimates and output the plots somewhere sensible. Find the means and standard deviations of your parameter estimates.
• Repeat the previous exercises with a mixture of two normal distributions as your model instead of a simple Gaussian (I've elected in my example code to only fit the x vector, but you can try and fit the full joint distribution if you want). Observe that it is a real challenge to get acceptable results in a reasonable length of time (this will be important in later chapters).
• Commit your work locally, then merge your development branch into main. Push main to your remote repository.

Notes: if you find yourself with a simulation that's running overlong, make sure to terminate it with Ctrl + C and then quit out of the interactive R session before closing the terminal. If you don't do this you can wind up with dangling Docker containers eating more and more CPU cores until your computer becomes unusable. If you do wind up in this state, sudo docker ps followed by sudo docker stop <container-id> will allow you to shut down any troublesome containers.

It can be difficult to tell whether you're making meaningful progress with some sets of hyperparameters. I thus highly recommend that you add a progress bar to the core rejection sampling scheme: the r-cli project has excellent support for them. These will, however, foul up any attempts to snapshot test the function, so you'll need a way to turn them off for testing purposes.

Finally, it's important that your tests run quickly: if they take too long to run, you'll get frustrated and won't run them. It's thus OK to generate small samples all the way through this: we're primarily interested in whether the software engineering is sound rather than statistical considerations, so you don't necessarily need to check whether or not the process converges to the correct value in your unit tests (that's a job for you to do as an analyst).

If you were able to complete all of this successfully, congratulations! You're now able to produce non-trivial, high-quality, well-engineered code in R. If you struggled a bit, read through the the example code base that I've provided carefully and try again: stepping through things bit by bit and not progressing until you're sure you understand the first thing is the best way to progress here. In the next chapter we're going to talk about databases and SQL: a subject that's relevant to statisticians, but of first-order importance for engineers.