Lab 02 - R Essentials and Functions

PHS 7045: Advanced Programming — Fall 2026

Learning goals

  • Warm up by writing loops and functions, and see why pkg::fun() is worth typing.
  • Turn a one-off simulation into a parameterized function.
  • Compare an R-level loop against a genuinely vectorized implementation, and measure the difference instead of guessing.
  • Use labels (named vectors, named columns) everywhere it makes the code clearer.
  • Take a first step toward HW 1.

The second half uses Bayesian statistics, but the emphasis is on the coding. We will walk through everything needed for the estimation.

Warm-up: loops, functions, and beepr

Install the package beepr and run beepr::beep(). beepr::beep(k) can play k = 1 to 11 sounds (see ?beep for the list).

install.packages("beepr")
beepr::beep(0)

A. Loops and functions

  1. Write a loop that plays each sound in turn. Use Sys.sleep() to pause 2 seconds between calls to beep.

  2. Modify the loop to pause a random duration of time. You can set your own parameters for what “a random duration of time” means.

B. How can beepr::beep() be useful? What does this say about R being a scripting language?

C. After library(beepr), you can call beep() directly. Why might you still write beepr::beep()?

D. Write a function that takes two inputs, a numeric vector and sound. Inside it, write a loop that calculates the cumulative sum one element at a time (do not use cumsum()), and play a beep at the end using the sound input. The output should be a two-column matrix with the original vector and the cumulative sum. Check your answer against cumsum().

Note that question D is designed to practice working with loops and building up a result. One-at-a-time calculations are discouraged whenever you can avoid them; in practice you would just use cumsum().

This warm-up is not submitted. We discuss B and C together before moving on.

Part 1: The urn problem

From The Art of R Programming:

Urn 1 contains ten blue marbles and eight yellow ones. In urn 2, the mixture is six blue and six yellow. We draw a marble at random from urn 1, transfer it to urn 2, and then draw a marble at random from urn 2. What is the probability that the second marble is blue?

This is easy to find analytically, but we will use simulation.

  1. Write a first version any way you like: a loop, replicate(), whatever comes to mind. Then turn it into a function with these arguments and defaults:

    urn_prob <- function(nreps = 100000, b1 = 10, y1 = 8, b2 = 6, y2 = 6) {
  # your code
}

  2. Check your function against the analytic answer, \(\frac{b_1}{b_1+y_1}\cdot\frac{b_2+1}{b_2+y_2+1} + \frac{y_1}{b_1+y_1}\cdot\frac{b_2}{b_2+y_2+1}\). Does it agree to within Monte Carlo error? How would you quantify “within Monte Carlo error”?

  3. Post your function here.

  4. We will benchmark the class’s solutions together and discuss how to speed them up:

    library(bench)
bench::mark(<solutions>, relative = TRUE, check = FALSE)

  5. Write a second version that removes the R-level loop entirely: draw all nreps replicates at once with rbinom(). Benchmark it against your first version. Where did the time actually go?

Submit: your original function, the more efficient version, the benchmark, and a short paragraph on what you would do differently next time.

Part 2: Extending the urn problem

The point of parameterizing the function is that new questions become cheap.

  1. Sweep the urn compositions. Use expand.grid() to build a grid of scenarios, varying b1 over 5:15 and holding the other urns fixed, then estimate the probability for each. Return a labeled data frame with one row per scenario. Use a functional (Map(), mapply(), or apply() over the rows) rather than a nested loop, and plot the result.

  2. Generalize to \(k\) urns. Extend the function so the marble is transferred through a chain of urns: draw from urn 1, transfer to urn 2, draw from urn 2, transfer to urn 3, and so on, reporting the probability that the final draw is blue. Represent the urns as a labeled structure (a matrix with a blue and yellow column and one row per urn, or a list of named vectors) so the code reads clearly.

    Note where vectorization stops being possible: each transfer depends on the previous draw, so the chain is inherently sequential. What can still be vectorized across replicates?

  3. Add a Monte Carlo standard error. Have your function return not just the estimate but the standard error, \(\sqrt{\hat p (1-\hat p) / \texttt{nreps}}\), as a named vector. How many replicates do you need for a standard error below 0.001?

Part 3: Estimating the probability one treatment arm is better than another

Create a simulated dataset that enrolls 40 participants (\(i = 1, \dots, 40\)) equally randomized to 4 arms (\(t = 0, 1, 2, 3\)). Generate outcomes supposing the probability of success under each arm, \(p_t\), is 0.35.

Begin by creating a matrix with 10 rows and 4 columns. Label the columns to indicate each treatment arm and generate the observed outcomes under each treatment arm. The total number of observations and successes under arm \(t\) is, respectively, \(n_t\) and \(y_t\).

Compare the probability of success under each experimental arm to control, i.e., Pr(\(p_t > p_0\)).

We will use a Bayesian framework to estimate the distribution of \(p_t\) under each arm. With Bernoulli outcomes and a Beta(\(\alpha_t = 0.35, \beta_t = 0.65\)) prior on the success rate, the posterior is Pr(\(p_t \mid y_t\)) \(\sim\) Beta(\(\alpha_t + y_t\), \(\beta_t + n_t - y_t\)).

To estimate Pr(\(p_t > p_0\)), take many (say 1000) random draws from each arm’s posterior: rbeta(n = 1000, 0.35 + y_t, 0.65 + n_t - y_t). Use vectorization to compare how often a randomly generated experimental success rate exceeds a randomly generated control success rate.

Declare the trial a “success” if the maximum of Pr(\(p_t > p_0\)) \(> \delta = 0.025\).

Throughout, use labels wherever appropriate. Submit this as your lab assignment.

Next steps

With minor modifications, this is the setup for the first of the two trial designs in HW 1. If you have extra time, start on HW 1.