library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)
library(plotly)      # interactive figures and ALL 3-D graphics
library(data.table)
library(bench)       # honest microbenchmarking

How this chapter uses graphics

Every two-dimensional figure is drawn with ggplot2 and rendered statically. Immediately after each one, the equivalent plot_ly() code appears in a chunk marked eval=FALSE, echo=TRUE.

Every three-dimensional figure is drawn with plot_ly() and evaluated. Computational performance is intrinsically multi-variable — speedup depends on both the parallel fraction and the core count, attainable throughput on both arithmetic intensity and hardware, sketch accuracy on both memory budget and stream length. Each is a surface, and the interaction is the whole point.

On portability. This chapter sets no working directory and contains no absolute paths. Every file is written under tempdir() and every network call is wrapped so that an unavailable service produces a stated fallback rather than a failed build.

1 Learning objectives

After completing this chapter you will be able to:

Explain the difference between row-major and columnar storage, and predict which operations each favours.
Benchmark CSV, RDS, Parquet, and Arrow on read time, file size, and column-pruning cost.
Push predicates down to a database or query engine, and measure the reduction in data transferred.
Perform larger-than-memory analytics with an out-of-core query engine.
State and apply reservoir sampling, the count-min sketch, and HyperLogLog, with their error guarantees.
Evaluate a streaming model prequentially and detect concept drift.
Locate a computation on the roofline and decide whether it is compute-bound or memory-bound.
Derive Amdahl’s and Gustafson’s laws and use them to predict the ceiling on parallel speedup.
Model parallel overhead and compute the core count beyond which parallelism loses.
Compare summation algorithms on both speed and numerical accuracy, and explain why a faster reduction can be a worse one.

Estimated time: 9–12 hours including exercises. Prerequisites: Chapter 1 (the R toolchain), Chapter 3 (floating point and conditioning, §3.7), and Chapter 8 (the adjusted Rand index, §8.4.2). This chapter is the engineering layer beneath all of them.

2 PART I: FORMATS AND INGESTION

Every previous chapter began with a tidy data frame. Producing one is often where most of the work lives, and the choices made there set hard limits on what is computationally possible afterwards.

3 How data is laid out in memory

A two-dimensional table has to be written into one-dimensional memory, and there are two ways to do it.

\[ \textbf{Row-major: }\ \underbrace{x_{11}\,x_{12}\,x_{13}}_{\text{row }1}\ \underbrace{x_{21}\,x_{22}\,x_{23}}_{\text{row }2}\ \cdots \qquad \textbf{Column-major: }\ \underbrace{x_{11}\,x_{21}}_{\text{col }1}\ \underbrace{x_{12}\,x_{22}}_{\text{col }2}\ \cdots \]

R matrices and data frames are column-major. C arrays and CSV files are row-major. That single fact drives a great deal of practical performance.

Common misconception: “the storage layout is an implementation detail.” It determines which operations are cheap. Summing a column in a column-major layout reads contiguous memory, so each cache line fetched delivers eight useful doubles. Summing a row strides across the array, and every access may pull a fresh 64-byte cache line to use 8 bytes of it — an eightfold waste of memory bandwidth.

The consequence for analytics is direct: most analytical queries touch few columns and many rows, which is exactly the access pattern columnar storage is built for. Most transactional queries touch one row and many columns, which is what row storage is built for.

set.seed(11)
n_l <- 4000
M <- matrix(rnorm(n_l * n_l), n_l, n_l)   # R stores this COLUMN-major

bench::mark(
  `sum a column (contiguous)` = sum(M[, 2000]),
  `sum a row (strided)`       = sum(M[2000, ]),
  check = FALSE, iterations = 40
)[, c("expression", "median", "mem_alloc")]

Same number of additions, same number of doubles, and a substantial timing difference, entirely attributable to memory access pattern.

4 Format conversion

The rio package reads and writes some forty formats behind three verbs: import(), export(), and convert(). It infers the format from the file extension.

library(rio)

# Everything is written under tempdir(); nothing touches the working directory
work_dir <- file.path(tempdir(), "ch10"); dir.create(work_dir, showWarnings = FALSE)

nof1 <- dspa_try(
  rio::import("https://umich.instructure.com/files/1760330/download?download_frd=1",
              format = "dta"),
  fallback = data.frame(ID = rep(1:10, each = 4), Day = rep(1:4, 10),
                        Tx = rep(0:1, 20), SelfEff = rnorm(40, 30, 5),
                        SelfEff25 = rnorm(40, 5, 2), WPSS = rnorm(40, 0, 1),
                        SocSuppt = rnorm(40, 5, 1), PMss = rnorm(40, 2, 1),
                        PMss3 = rnorm(40, -1, 1), PhyAct = rnorm(40, 200, 60)),
  label = "N-of-1 Stata file")
str(nof1)

#> 'data.frame':    900 obs. of  10 variables:
#>  $ id       : num  1 1 1 1 1 1 1 1 1 1 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ day      : num  1 2 3 4 5 6 7 8 9 10 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ tx       : num  1 1 0 0 1 1 0 0 1 1 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ selfeff  : num  33 33 33 33 33 33 33 33 33 33 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ selfeff25: num  8 8 8 8 8 8 8 8 8 8 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ wpss     : num  0.97 -0.17 0.81 -0.41 0.59 -1.16 0.3 -0.34 -0.74 -0.38 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ scssuppt : num  5 3.87 4.84 3.62 4.62 2.87 4.33 3.69 3.29 3.66 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ pmss     : num  4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 4.03 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ pmss3    : num  1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 1.03 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"
#>  $ qhyact   : num  53 73 23 36 21 0 21 0 73 114 ...
#>   ..- attr(*, "format.stata")= chr "%12.0g"

paths <- file.path(work_dir, c("nof1.csv", "nof1.rds", "nof1.xlsx", "nof1.parquet"))
names(paths) <- c("csv", "rds", "xlsx", "parquet")

rio::export(nof1, paths[["csv"]])
rio::export(nof1, paths[["rds"]])
ok_xlsx <- !is.null(dspa_try(rio::export(nof1, paths[["xlsx"]]),
                             label = "xlsx writer"))

# convert() is import() + export() in one call
rio::convert(paths[["csv"]], file.path(work_dir, "nof1_converted.tsv"))

data.frame(file = basename(list.files(work_dir)),
           bytes = file.size(list.files(work_dir, full.names = TRUE)))

5 Columnar formats

CSV is universal and slow: it is text, so every number must be parsed; it is row-major, so reading one column requires scanning every byte; and it carries no types, so they must be inferred.

Apache Parquet stores data by column, with per-column compression and per-chunk statistics (min, max, null count). Apache Arrow specifies an in-memory columnar layout that many languages can read without copying. Feather is Arrow written to disk.

Three properties follow, and they are what make columnar formats the modern default:

Property	Consequence
Column pruning	Reading 3 of 200 columns costs \(\approx 3/200\) of the I/O
Predicate pushdown	Row groups whose min/max exclude the filter are skipped unread
Type-aware compression	Run-length, dictionary, and delta encodings beat generic gzip on columnar data

library(arrow)
set.seed(21)
n_c <- 2e5
big <- data.frame(
  id       = seq_len(n_c),
  group    = sample(letters[1:8], n_c, TRUE),          # low cardinality: dictionary-encodes
  ts       = as.Date("2020-01-01") + sample(0:1500, n_c, TRUE),
  x        = rnorm(n_c), y = rnorm(n_c), z = rnorm(n_c),
  flag     = sample(c(TRUE, FALSE), n_c, TRUE)
)
for (j in 1:20) big[[paste0("v", j)]] <- rnorm(n_c)     # 27 columns total

f_csv <- file.path(tempdir(), "big.csv")
f_rds <- file.path(tempdir(), "big.rds")
f_pq  <- file.path(tempdir(), "big.parquet")
f_fth <- file.path(tempdir(), "big.feather")

data.table::fwrite(big, f_csv)
saveRDS(big, f_rds, compress = TRUE)
arrow::write_parquet(big, f_pq)
arrow::write_feather(big, f_fth)

sizes <- data.frame(
  format = c("CSV", "RDS (gzip)", "Parquet", "Feather"),
  MB = round(file.size(c(f_csv, f_rds, f_pq, f_fth)) / 1e6, 2))
sizes$ratio_vs_csv <- round(sizes$MB / sizes$MB[1], 3)
sizes

b_full <- bench::mark(
  CSV     = data.table::fread(f_csv, showProgress = FALSE),
  RDS     = readRDS(f_rds),
  Parquet = arrow::read_parquet(f_pq),
  Feather = arrow::read_feather(f_fth),
  check = FALSE, iterations = 5, memory = FALSE)

b_cols <- bench::mark(
  `CSV, 3 cols`     = data.table::fread(f_csv, select = c("id", "group", "x"),
                                        showProgress = FALSE),
  `Parquet, 3 cols` = arrow::read_parquet(f_pq, col_select = c("id", "group", "x")),
  `Feather, 3 cols` = arrow::read_feather(f_fth, col_select = c("id", "group", "x")),
  check = FALSE, iterations = 5, memory = FALSE)

rbind(
  data.frame(task = "read all 27 columns", format = as.character(b_full$expression),
             median_sec = round(as.numeric(b_full$median), 4)),
  data.frame(task = "read 3 of 27 columns", format = as.character(b_cols$expression),
             median_sec = round(as.numeric(b_cols$median), 4)))

Read the second block against the first. For CSV, selecting three columns saves little, the parser must still scan every byte of every row to find the delimiters. For Parquet and Feather the saving is close to proportional, because the columns are stored separately and the unread ones are never touched.

bind_rows(
  data.frame(format = as.character(b_full$expression),
             seconds = as.numeric(b_full$median), task = "All 27 columns"),
  data.frame(format = sub(", 3 cols", "", as.character(b_cols$expression)),
             seconds = as.numeric(b_cols$median), task = "3 of 27 columns")) |>
  ggplot(aes(reorder(format, seconds), seconds, fill = task)) +
  geom_col(position = "dodge", width = 0.7) +
  coord_flip() +
  scale_fill_manual(values = c("All 27 columns" = "#9EB4C8",
                               "3 of 27 columns" = "#3B7DD8")) +
  labs(title = "Column pruning is free in columnar formats and expensive in CSV",
       subtitle = "Same data, same columns requested. Only the storage layout differs",
       x = NULL, y = "Median read time (seconds)", fill = NULL) +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(x = as.character(b_full$expression), y = as.numeric(b_full$median),
        type = "bar", name = "All columns") |>
  add_trace(x = sub(", 3 cols", "", as.character(b_cols$expression)),
            y = as.numeric(b_cols$median), name = "3 columns") |>
  layout(title = "Read time by format and column selection",
         xaxis = list(title = ""), yaxis = list(title = "Seconds"),
         barmode = "group")

The advantage is a surface over (rows, columns):

row_grid <- c(2e4, 5e4, 1e5, 2e5)
col_grid <- c(5, 10, 20, 40)

read_time <- function(nr, nc, fmt) {
  set.seed(nr + nc)
  d <- as.data.frame(matrix(rnorm(nr * nc), nr, nc))
  f <- file.path(tempdir(), sprintf("t_%d_%d.%s", nr, nc, fmt))
  if (fmt == "csv") data.table::fwrite(d, f) else arrow::write_parquet(d, f)
  sel <- names(d)[1:min(3, nc)]
  t0 <- Sys.time()
  if (fmt == "csv") data.table::fread(f, select = sel, showProgress = FALSE)
  else arrow::read_parquet(f, col_select = all_of(sel))
  as.numeric(difftime(Sys.time(), t0, units = "secs"))
}

Zcsv <- outer(col_grid, row_grid, Vectorize(\(c_, r_) read_time(r_, c_, "csv")))
Zpq  <- outer(col_grid, row_grid, Vectorize(\(c_, r_) read_time(r_, c_, "parquet")))

plot_ly() |>
  add_surface(x = row_grid, y = col_grid, z = Zcsv, opacity = 0.9,
              showscale = FALSE, colorscale = "Reds", name = "CSV") |>
  add_surface(x = row_grid, y = col_grid, z = Zpq, opacity = 0.9,
              showscale = FALSE, colorscale = "Blues", name = "Parquet") |>
  layout(title = "Time to read 3 columns: CSV (red, upper) vs. Parquet (blue, lower)",
         scene = list(xaxis = list(title = "Rows"),
                      yaxis = list(title = "Total columns"),
                      zaxis = list(title = "Seconds")))

Rotate along the columns axis. The CSV surface climbs steeply as total column count grows, the parser scans every one regardless of what was requested, while the Parquet surface stays nearly flat, because the unread columns are never touched.

6 Relational databases and predicate pushdown

A query is an instruction sent to a database describing what is wanted; the engine decides how to obtain it. That separation is the entire value proposition, and the single most important consequence is predicate pushdown: filtering and aggregating happen where the data lives, so only the result crosses the wire.

# BRFSS_MODE = "sample" builds a synthetic survey table with the same shape as
# the CDC BRFSS extract; "full" reproduces the 315 MB ingest. Neither writes
# outside tempdir().
make_survey <- function(n = 50000, seed = 31) {
  set.seed(seed)
  race_lv <- c("White", "Black", "Asian", "AmIndian", "Hawaiian",
               "Other", "Multiracial", "Hispanic", "Refused")
  d <- data.frame(
    year     = sample(c(2013L, 2015L), n, TRUE),
    x_race   = factor(sample(race_lv, n, TRUE,
                             prob = c(.62,.12,.05,.02,.01,.03,.03,.11,.01)),
                      levels = race_lv),
    hlthpln1 = sample(1:2, n, TRUE, prob = c(0.86, 0.14)),
    poorhlth = pmax(0, round(rnorm(n, 4, 6))),
    genhlth  = sample(1:5, n, TRUE, prob = c(.18,.32,.30,.15,.05)),
    age      = sample(18:95, n, TRUE))
  for (j in 1:24) d[[sprintf("q%02d", j)]] <- sample(1:9, n, TRUE)
  d
}

survey <- if (BRFSS_MODE == "full" && ONLINE) {
  zp <- dspa_download("https://www.socr.umich.edu/data/DSPA/BRFSS_2013_2014_2015.zip",
                      "BRFSS.zip", timeout = 900)
  dspa_try({
    ex <- file.path(tempdir(), "brfss"); dir.create(ex, showWarnings = FALSE)
    members <- unzip(zp, list = TRUE)$Name          # select BY NAME, not position
    keep <- grep("2013|2015", members, value = TRUE)
    unzip(zp, files = keep, exdir = ex)             # extract ONCE, into tempdir
    on.exit(unlink(ex, recursive = TRUE), add = TRUE)
    do.call(rbind, lapply(keep, \(f) Hmisc::sasxport.get(file.path(ex, f))))
  }, fallback = make_survey(), label = "BRFSS archive")
} else make_survey()

c(mode = BRFSS_MODE, rows = nrow(survey), columns = ncol(survey),
  MB_in_memory = round(as.numeric(object.size(survey)) / 1e6, 1))

#>         mode         rows      columns MB_in_memory 
#>     "sample"      "50000"         "30"        "6.2"

# The outcome is created BEFORE it is summarized
survey$has_plan <- survey$hlthpln1 == 1
c(coverage_rate = round(mean(survey$has_plan), 4))

#> coverage_rate 
#>        0.8604

summary(survey$has_plan)

#>    Mode   FALSE    TRUE 
#> logical    6981   43019

round(prop.table(table(survey$x_race)), 4)

#> 
#>       White       Black       Asian    AmIndian    Hawaiian       Other 
#>      0.6264      0.1165      0.0494      0.0195      0.0103      0.0292 
#> Multiracial    Hispanic     Refused 
#>      0.0298      0.1091      0.0098

t_glm <- system.time(
  fit_hcp <- glm(has_plan ~ x_race, data = survey, family = binomial()))
c(seconds = round(t_glm[["elapsed"]], 2))

#> seconds 
#>    0.09

# Coefficients are LOG-odds ratios relative to the reference level.
# Exponentiating gives odds ratios (Chapter 9, Section 9.8).
or_tab <- data.frame(
  level = names(coef(fit_hcp))[-1],
  log_OR = round(coef(fit_hcp)[-1], 4),
  OR = round(exp(coef(fit_hcp)[-1]), 4),
  CI_low = round(exp(confint.default(fit_hcp)[-1, 1]), 4),
  CI_high = round(exp(confint.default(fit_hcp)[-1, 2]), 4))
head(or_tab, 5)

6.1 Doing the work in the database

library(DBI)
library(duckdb)

con <- dbConnect(duckdb::duckdb(), dbdir = ":memory:")

dbWriteTable(con, "survey", survey)
dbWriteTable(con, "USArrests", tibble::rownames_to_column(USArrests, "state"))
dbListTables(con)

#> [1] "USArrests" "survey"

Common misconception: “the point of a database is that it stores the data.” The point is that it computes. SELECT * FROM t transfers the whole table into R and throws away every advantage the engine offers, the database becomes an expensive file.

What a query engine gives you is predicate pushdown (filter before transfer), projection pushdown (send only the requested columns), aggregation pushdown (compute the summary server-side and return a few rows), and an optimizer that reorders operations. Push the computation down and the data that crosses the wire is the answer, not the input.

# Aggregation done IN the database: the result is a handful of rows
  t_push <- system.time(
    agg_db <- dbGetQuery(con, "
      SELECT year, hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
      FROM survey
      WHERE age BETWEEN 45 AND 64
      GROUP BY year, hlthpln1
      ORDER BY year, hlthpln1"))
  
  # Everything pulled into R first, then aggregated
  t_pull <- system.time({
    raw <- dbGetQuery(con, "SELECT * FROM survey")
    agg_r <- raw |> filter(age >= 45, age <= 64) |>
      summarise(n = dplyr::n(), mean_poorhlth = mean(poorhlth),
                .by = c(year, hlthpln1)) |> arrange(year, hlthpln1)
  })
  
  agg_db

  c(pushdown_seconds = round(t_push[["elapsed"]], 3),
    pull_then_aggregate_seconds = round(t_pull[["elapsed"]], 3),
    rows_transferred_pushdown = nrow(agg_db),
    rows_transferred_pull = nrow(raw),
    transfer_ratio = round(nrow(raw) / nrow(agg_db)))

#>            pushdown_seconds pull_then_aggregate_seconds 
#>                           0                           0 
#>   rows_transferred_pushdown       rows_transferred_pull 
#>                           4                       50000 
#>              transfer_ratio 
#>                       12500

The pushdown query transfers a handful of rows; the pull-then-aggregate version transfers the entire table. The ratio is the saving, and it grows linearly with the data.

library(dbplyr)

# dplyr verbs against a database table are TRANSLATED to SQL and executed there.
# Nothing is materialized until collect().
tbl_survey <- tbl(con, "survey")

q <- tbl_survey |>
  filter(age >= 45, age <= 64) |>
  summarise(n = n(), mean_poorhlth = mean(poorhlth, na.rm = TRUE),
            .by = c(year, hlthpln1))

# Inspect the generated SQL before running anything
show_query(q)

#> <SQL>
#> SELECT "year", hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
#> FROM (
#>   SELECT survey.*
#>   FROM survey
#>   WHERE (age >= 45.0) AND (age <= 64.0)
#> ) q01
#> GROUP BY "year", hlthpln1

collect(q) |> arrange(year, hlthpln1)

show_query() is the habit worth forming: it makes visible exactly what was pushed down and what was not.

6.2 Out-of-core analytics

When a dataset exceeds RAM, the classical answers were disk-backed data structures. The modern answer is a query engine that streams over columnar files, keeping only the working set in memory.

# Write a partitioned Parquet dataset -- larger than we want to hold at once
ds_dir <- file.path(tempdir(), "survey_ds")
arrow::write_dataset(survey, ds_dir, partitioning = "year", format = "parquet")

list.files(ds_dir, recursive = TRUE)[1:4]

#> [1] "year=2013/part-0.parquet" "year=2015/part-0.parquet"
#> [3] NA                         NA

c(total_MB_on_disk = round(sum(file.size(list.files(ds_dir, recursive = TRUE,
                                                    full.names = TRUE))) / 1e6, 2))

#> total_MB_on_disk 
#>             0.77

# Query the dataset WITHOUT loading it: arrow streams over the files
ds <- arrow::open_dataset(ds_dir)
class(ds)

#> [1] "FileSystemDataset" "Dataset"           "ArrowObject"      
#> [4] "R6"

res_arrow <- ds |>
  filter(year == 2015, age >= 45, age <= 64) |>
  summarise(n = n(), mean_poorhlth = mean(poorhlth), .by = hlthpln1) |>
  collect()
res_arrow

# duckdb can query Parquet files directly, with partition pruning
sql <- sprintf("
  SELECT hlthpln1, COUNT(*) AS n, AVG(poorhlth) AS mean_poorhlth
  FROM read_parquet('%s/**/*.parquet', hive_partitioning = true)
  WHERE year = 2015 AND age BETWEEN 45 AND 64
  GROUP BY hlthpln1", ds_dir)
dbGetQuery(con, sql)

c(note = "The 2013 partition is never opened: Hive partitioning prunes it by directory")

#>                                                                           note 
#> "The 2013 partition is never opened: Hive partitioning prunes it by directory"

Partition pruning happens before any file is read. Because year is encoded in the directory structure, a filter on year eliminates whole directories, and a filter on age then uses Parquet’s per-row-group min/max statistics to skip chunks. Neither requires decompressing the data.

# A dataset deliberately larger than the R session should hold at once
big_dir <- file.path(tempdir(), "big_ds"); dir.create(big_dir, showWarnings = FALSE)
set.seed(41)
for (part in 1:6) {
  d <- data.frame(part = part,
                  g = sample(letters[1:10], 3e5, TRUE),
                  v = rnorm(3e5), w = rnorm(3e5))
  arrow::write_parquet(d, file.path(big_dir, sprintf("part-%d.parquet", part)))
}
c(rows_total = 6 * 3e5,
  MB_on_disk = round(sum(file.size(list.files(big_dir, full.names = TRUE))) / 1e6, 1))

#> rows_total MB_on_disk 
#>    1800000         33

t_ooc <- system.time(
  out <- arrow::open_dataset(big_dir) |>
    filter(v > 1) |>
    summarise(n = n(), mean_w = mean(w), .by = g) |>
    collect() |> arrange(g))
out

c(seconds = round(t_ooc[["elapsed"]], 3),
  peak_R_objects_MB = round(as.numeric(object.size(out)) / 1e6, 4))

#>           seconds peak_R_objects_MB 
#>             0.060             0.002

The result occupies kilobytes; the source occupies tens of megabytes on disk and was never fully materialized. This is the pattern that replaces disk-backed data frames: leave the data in a columnar file, express the query in dplyr, and let the engine stream.

7 Web data

7.1 JSON

library(jsonlite)
nof1_json <- dspa_try(
  jsonlite::fromJSON("https://umich.instructure.com/files/1760327/download?download_frd=1"),
  fallback = nof1, label = "N-of-1 JSON")
class(nof1_json); dim(nof1_json)

#> [1] "data.frame"

#> [1] 900  10

head(nof1_json, 3)

JSON’s {key: value} nesting maps naturally onto R lists; jsonlite::fromJSON(simplifyDataFrame = TRUE) flattens a regular array of objects into a data frame. When the structure is irregular, keep it as a list and use purrr::map_* to extract fields.

7.2 XML and scraping

library(rvest); library(xml2)

socr <- dspa_try(rvest::read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data"),
                 label = "SOCR wiki")
if (!is.null(socr)) {
  cat("Title:", socr |> html_element("head title") |> html_text(), "\n")
  paras <- socr |> html_elements("p") |> html_text()
  cat("Paragraphs found:", length(paras), "\n")
  cat(substr(paste(paras, collapse = " "), 1, 220), "...\n")
}

#> Title: SOCR Data - SOCR 
#> Paragraphs found: 29 
#> The links below contain a number of datasets that may be used for demonstration purposes in probability and statistics education. There are two types of data - simulated (computer-generated using random sampling) and obs ...

Check the terms of service and robots.txt before scraping, and rate-limit requests. polite::bow() handles both automatically. Scraping is a last resort: if an API exists, use it, because it is stable and the HTML is not.

7.3 Semantic queries with SPARQL

RDF represents knowledge as (subject, predicate, object) triples, and SPARQL queries patterns over those triples. Wikidata exposes a public endpoint.

library(WikidataQueryServiceR)

tb_query <- '
PREFIX wd:        <http://www.wikidata.org/entity/>
PREFIX wdt:       <http://www.wikidata.org/prop/direct/>
PREFIX rdfs:      <http://www.w3.org/2000/01/rdf-schema#>
PREFIX p:         <http://www.wikidata.org/prop/>
PREFIX qualifier: <http://www.wikidata.org/prop/qualifier/>
PREFIX statement: <http://www.wikidata.org/prop/statement/>

SELECT DISTINCT ?countryLabel ?ISO3Code ?latlon ?prevalence ?year
WHERE {
  wd:Q12204 p:P1193 ?prevalenceStatement .          # Q12204 = tuberculosis
  ?prevalenceStatement qualifier:P17  ?country ;
                       qualifier:P585 ?year ;
                       statement:P1193 ?prevalence .
  ?country wdt:P625  ?latlon ;
           rdfs:label ?countryLabel ;
           wdt:P298  ?ISO3Code .
  FILTER (lang(?countryLabel) = "en")
}
ORDER BY DESC(?prevalence)                           # ORDER BY a BOUND variable
'

tb <- dspa_try(WikidataQueryServiceR::query_wikidata(tb_query),
               label = "Wikidata SPARQL endpoint")
if (!is.null(tb)) { cat("Rows returned:", nrow(tb), "\n"); head(tb, 4) }

#> Rows returned: 6

The ORDER BY clause sorts on ?prevalence, which the WHERE block binds. Ordering by an unbound variable is a silent no-op, SPARQL treats unbound values as mutually incomparable, so every sort key must appear in the pattern.

if (!is.null(tb) && nrow(tb) > 0) {
  # WKT is "Point(longitude latitude)". Split, strip the affixes, and CONVERT
  # to numeric -- plot_geo requires numbers, not the character strings that
  # gsub leaves behind.
  geo <- tb |>
    tidyr::separate(latlon, into = c("long", "lat"), sep = " ", remove = FALSE) |>
    mutate(long = as.numeric(gsub("Point\\(", "", long)),
           lat  = as.numeric(gsub("\\)", "", lat)),
           prevalence = as.numeric(prevalence)) |>
    filter(is.finite(long), is.finite(lat), is.finite(prevalence))

  # Quantile binning is guarded: Wikidata coverage is sparse and may return
  # too few distinct values for four bins.
  n_distinct_prev <- length(unique(geo$prevalence))
  n_bins <- max(2, min(4, n_distinct_prev - 1))
  brk <- unique(quantile(geo$prevalence, seq(0, 1, length.out = n_bins + 1)))
  geo$band <- if (length(brk) >= 3)
    cut(geo$prevalence, breaks = brk, include.lowest = TRUE,
        labels = paste0("Q", seq_len(length(brk) - 1))) else factor("all")

  c(countries = nrow(geo), distinct_prevalences = n_distinct_prev, bins = n_bins)
  head(geo[, c("countryLabel", "ISO3Code", "long", "lat", "prevalence", "band")], 5)
}

if (exists("geo") && !is.null(geo) && nrow(geo) > 0) {
  ge <- list(scope = "world", showland = TRUE, landcolor = toRGB("gray92"),
             countrycolor = toRGB("white"), countrywidth = 0.5)
  plot_geo(geo, lon = ~long, lat = ~lat, color = ~band,
           text = ~paste0(countryLabel, " (", ISO3Code, "): ", prevalence),
           marker = list(size = 12, line = list(width = 0.5))) |>
    layout(geo = ge, title = "Tuberculosis prevalence reported in Wikidata")
}

library(maps)
cities <- maps::world.cities |>
  filter(pop > 5e5) |>
  mutate(label = sprintf("%s, %s: %.1fM", name, country.etc, pop / 1e6),
         # Marker AREA proportional to population, so radius scales as sqrt()
         size = sqrt(pop / max(pop)) * 22 + 2)

ge <- list(scope = "world", showland = TRUE, landcolor = toRGB("gray92"),
           countrycolor = toRGB("white"), countrywidth = 0.5)

plot_geo(cities, lon = ~long, lat = ~lat, text = ~label,
         marker = ~list(size = size, opacity = 0.6,
                        line = list(width = 0.2))) |>
  layout(geo = ge, title = "World cities above 500,000 population")

Marker area is proportional to population, so the radius scales as \(\sqrt{\text{pop}}\). Scaling the radius linearly with the value would exaggerate large cities quadratically, the same perceptual point made in Chapter 2, §2.10.

7.4 Random number services

rng_url <- paste0("https://www.random.org/integers/",
                  "?num=300&min=100&max=200&col=3&base=10&format=plain&rnd=new")
true_rng <- dspa_try(utils::read.table(rng_url),
                     fallback = as.data.frame(matrix(sample(100:200, 300, TRUE), ncol = 3)),
                     label = "random.org")
head(true_rng, 3)

Common misconception: “true random numbers are better for simulation.” For Monte Carlo work they are usually worse, because they are not reproducible. The property that matters is not physical randomness but statistical adequacy plus a seed: a modern PRNG (Mersenne Twister, PCG64) passes the standard test batteries and regenerates the identical stream from set.seed(), which is what makes a simulation checkable.

Physical entropy matters for cryptography, where predictability is the threat model, and for seeding a PRNG. It is not an upgrade for a simulation study.

8 PART II: NETWORK DATA

9 Graph representations and their cost

A graph \(G=(V,E)\) with \(|V|=n\) vertices and \(|E|=m\) edges admits three standard representations, and the choice governs everything downstream.

Representation	Memory	Edge query \((u,v)\)	Enumerate neighbours of \(u\)
Adjacency matrix	\(\Theta(n^2)\)	\(O(1)\)	\(O(n)\)
Adjacency list	\(\Theta(n+m)\)	\(O(\deg u)\)	\(O(\deg u)\)
Edge list	\(\Theta(m)\)	\(O(m)\)	\(O(m)\)

Real networks are sparse: \(m=O(n)\) or \(O(n\log n)\), far below the \(\binom n2\) maximum. For the Facebook ego network below, \(n\approx4{,}039\) and \(m\approx88{,}000\), so the adjacency matrix would hold \(1.6\times10^7\) entries of which 99.5% are zero.

graph_memory <- function(n, m) {
  c(n = n, m = m,
    density = signif(2 * m / (n * (n - 1)), 3),
    matrix_MB = round(8 * n^2 / 1e6, 1),
    list_MB   = round(8 * (n + 2 * m) / 1e6, 3),
    ratio = round((8 * n^2) / (8 * (n + 2 * m)), 1))
}
as.data.frame(rbind(
  `Facebook ego net` = graph_memory(4039, 88234),
  `Les Miserables`   = graph_memory(77, 254),
  `Hypothetical 1e5` = graph_memory(1e5, 1e6)))

At \(n=10^5\) the dense matrix needs 80 GB and the adjacency list needs 17 MB. Sparse representation is not an optimization; it is the difference between possible and impossible.

library(igraph)
g_toy <- igraph::make_graph(c(1,2, 1,3, 2,3, 3,4), n = 10, directed = TRUE)
plot(g_toy, vertex.size = 22, vertex.color = "#9EC5E8",
     edge.arrow.size = 0.5, main = "Four directed edges among ten vertices")

c(vertices = gorder(g_toy), edges = gsize(g_toy),
  isolated = sum(degree(g_toy, mode = "all") == 0))

#> vertices    edges isolated 
#>       10        4        6

Vertices 5–10 appear because n = 10 was declared; they have degree zero because no edge mentions them. Isolated vertices are a common artifact of declaring a vertex count larger than the edge list requires.

10 A real network

fb_path <- dspa_download(
  "https://umich.instructure.com/files/2854431/download?download_frd=1",
  "facebook_edges.txt")

edges <- dspa_try({
  e <- data.table::fread(fb_path, header = FALSE, col.names = c("from", "to"))
  as.matrix(e)
}, fallback = {
  set.seed(51)
  igraph::as_edgelist(igraph::sample_pa(1500, m = 4, directed = FALSE))
}, label = "Facebook ego network")

# The degenerate vertex is IDENTIFIED, not removed by a magic row count
c(min_id = min(edges), rows_touching_min = sum(edges[, 1] == min(edges) |
                                                edges[, 2] == min(edges)))

#>            min_id rows_touching_min 
#>                 0               347

edges <- edges[edges[, 1] != min(edges) & edges[, 2] != min(edges), , drop = FALSE]

g_fb <- igraph::graph_from_edgelist(edges, directed = FALSE)
c(vertices = gorder(g_fb), edges = gsize(g_fb),
  density = signif(edge_density(g_fb), 4),
  memory_MB = round(as.numeric(object.size(g_fb)) / 1e6, 2))

#>   vertices      edges    density  memory_MB 
#> 4.0380e+03 8.7887e+04 1.0780e-02 1.4100e+00

10.1 Centrality and what it costs

Degree centrality counts a vertex’s edges, \(O(1)\) per vertex from an adjacency list.

Betweenness centrality measures how often a vertex lies on shortest paths:

\[g(v)=\sum_{s\ne v\ne t}\frac{\sigma_{st}(v)}{\sigma_{st}},\]

with \(\sigma_{st}\) the number of shortest \(s\)–\(t\) paths and \(\sigma_{st}(v)\) the number passing through \(v\).

Betweenness is expensive. The naive definition suggests enumerating all \(\binom n2\) pairs. Brandes’ algorithm (2001) reduces this to \(O(nm)\) for unweighted graphs and \(O(nm+n^2\log n)\) for weighted ones, a genuine breakthrough, and still \(O(nm)\). For the Facebook network that is \(4{,}038\times87{,}887\approx3.5\times10^8\) operations.

For larger graphs, sample the source vertices: estimating from \(k\ll n\) randomly chosen sources gives \(O(km)\) with an error that shrinks as \(O(1/\sqrt k)\).

deg <- igraph::degree(g_fb)
hub <- which.max(deg)
c(n_vertices = length(deg),
  max_degree = max(deg),
  hub_vertex = as.integer(names(hub)),
  mean_degree = round(mean(deg), 2),
  median_degree = median(deg))

#>    n_vertices    max_degree   mean_degree median_degree 
#>       4038.00       1044.00         43.53         25.00

# Exact betweenness on a graph this size is O(nm); sample sources instead
t_exact <- system.time(bw_exact <- igraph::betweenness(g_fb, normalized = FALSE))
set.seed(53)
t_approx <- system.time(
  bw_approx <- igraph::betweenness(g_fb, normalized = FALSE,
                                   cutoff = 4))     # bounded-length paths

c(exact_seconds = round(t_exact[["elapsed"]], 2),
  cutoff_seconds = round(t_approx[["elapsed"]], 2),
  correlation = round(cor(bw_exact, bw_approx, method = "spearman"), 4),
  top_vertex_exact = as.integer(names(which.max(bw_exact))),
  top_vertex_cutoff = as.integer(names(which.max(bw_approx))))

#>  exact_seconds cutoff_seconds    correlation 
#>         1.3600         1.0200         0.9829

The which.max calls compute the maximum over all vertices. Reading a maximum off a ten-element window would establish nothing.

dd <- data.frame(degree = as.integer(names(table(deg))),
                 count = as.integer(table(deg))) |> filter(degree > 0)

ggplot(dd, aes(degree, count)) +
  geom_point(alpha = 0.6, size = 1.3, color = "steelblue") +
  scale_x_log10() + scale_y_log10() +
  labs(title = "Degree distribution of the ego network",
       subtitle = "Approximately linear on log-log axes: a heavy tail, with a few very high-degree hubs",
       x = "Degree (log)", y = "Number of vertices (log)") +
  theme_dspa()

# --- Interactive equivalent ------------------------------------------------
plot_ly(dd, x = ~degree, y = ~count, type = "scatter", mode = "markers") |>
  layout(title = "Degree distribution",
         xaxis = list(title = "Degree", type = "log"),
         yaxis = list(title = "Count", type = "log"))

10.2 Community detection

set.seed(55)
t_louvain <- system.time(comm <- igraph::cluster_louvain(g_fb))
c(communities = length(comm),
  modularity = round(modularity(comm), 4),
  seconds = round(t_louvain[["elapsed"]], 2),
  largest = max(sizes(comm)), smallest = min(sizes(comm)))

#> communities  modularity     seconds     largest    smallest 
#>     35.0000      0.8341      0.0300    548.0000      1.0000

head(sort(sizes(comm), decreasing = TRUE), 8)

#> Community sizes
#>  28  30   9  14  16  25   1  32 
#> 548 535 459 446 423 323 238 237

Modularity compares the within-community edge fraction against what a degree-preserving random graph would produce:

\[Q=\frac{1}{2m}\sum_{i,j}\left(A_{ij}-\frac{k_ik_j}{2m}\right)\delta(c_i,c_j) \ \in\ \left[-\tfrac12,\,1\right).\]

The subtracted term is the configuration-model null, the same null-comparison logic as the gap statistic (Chapter 8, §8.7.3) and the permutation test for association rules (Chapter 7, §7.18.3). Louvain is \(O(m\log n)\) in practice and greedy, so it finds a local optimum; Leiden (cluster_leiden) repairs a known Louvain defect in which communities can be internally disconnected.

sz <- as.integer(sizes(comm))
ggplot(data.frame(size = sz), aes(size)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  scale_x_log10() +
  labs(title = sprintf("%d communities detected (modularity %.3f)",
                       length(comm), modularity(comm)),
       subtitle = "Sizes span orders of magnitude -- a few large communities and many small ones",
       x = "Community size (log scale)", y = "Count") +
  theme_dspa()

11 PART III: DATA STREAMS

12 The streaming model

A data stream is an ordered, potentially unbounded sequence

\[Y=\{y_1,y_2,\dots,y_t,\dots\}\]

arriving in order and too large to store. High-frequency examples include ICU EEG telemetry, market data, and surface weather observations.

The model imposes three constraints that reshape every algorithm:

One pass. Each item is seen once; there is no rewinding.
Sublinear memory. Typically \(O(\mathrm{polylog}\,n)\), never \(O(n)\).
Anytime answers. A current estimate must be available at any moment.

Common misconception: “streaming just means processing data in chunks.” Chunked batch processing still allows multiple passes and \(O(n)\) memory. Streaming forbids both, and that forces genuinely different algorithms: exact answers become impossible and are replaced by sketches, compact summaries with provable \((\varepsilon,\delta)\) error guarantees.

The trade is explicit and quantified: you accept a small, bounded error in exchange for memory that does not grow with the stream.

13 Sketches

13.1 Reservoir sampling

Problem. Draw a uniform sample of size \(k\) from a stream of unknown length, in one pass and \(O(k)\) memory.

Algorithm R (Vitter, 1985): keep the first \(k\) items. For \(t>k\), generate \(j\sim\mathrm{Unif}\{1,\dots,t\}\) and if \(j\le k\), replace reservoir slot \(j\) with item \(t\).

Theorem. After processing \(n\) items, every item is in the reservoir with probability exactly \(k/n\).

Proof sketch. By induction. Item \(t>k\) enters with probability \(k/t\). An item already present survives step \(t\) with probability \(1-\frac{k}{t}\cdot\frac1k=\frac{t-1}{t}\). Telescoping from \(t+1\) to \(n\) gives survival probability \(\frac{t}{n}\), so the total is \(\frac kt\cdot\frac tn=\frac kn\). \(\blacksquare\)

reservoir_sample <- function(stream, k) {
  res <- integer(k)
  for (t in seq_along(stream)) {
    if (t <= k) res[t] <- stream[t]
    else { j <- sample.int(t, 1); if (j <= k) res[j] <- stream[t] }
  }
  res
}

set.seed(61)
n_r <- 2000; k_r <- 50
counts <- tabulate(unlist(replicate(3000, reservoir_sample(1:n_r, k_r),
                                    simplify = FALSE)), nbins = n_r)
emp <- counts / 3000

c(theoretical_inclusion_prob = k_r / n_r,
  mean_empirical = round(mean(emp), 5),
  sd_empirical = round(sd(emp), 5),
  memory = sprintf("O(k) = %d items, regardless of stream length", k_r))

#>                     theoretical_inclusion_prob 
#>                                        "0.025" 
#>                                 mean_empirical 
#>                                        "0.025" 
#>                                   sd_empirical 
#>                                      "0.00281" 
#>                                         memory 
#> "O(k) = 50 items, regardless of stream length"

ggplot(data.frame(position = seq_len(n_r), p = emp), aes(position, p)) +
  geom_point(alpha = 0.25, size = 0.6, color = "steelblue") +
  geom_hline(yintercept = k_r / n_r, color = "firebrick", linewidth = 1) +
  labs(title = "Reservoir sampling gives every item the same inclusion probability",
       subtitle = sprintf("Red line: theoretical k/n = %.3f. Position along the stream has no effect",
                          k_r / n_r),
       x = "Position in the stream", y = "Empirical inclusion probability") +
  theme_dspa()

The empirical probabilities scatter around \(k/n\) with no trend in position — early and late items are equally likely, which is exactly what uniformity requires and what a naive “keep the first \(k\)” or “keep the last \(k\)” scheme would violate.

13.2 The count-min sketch

Problem. Estimate the frequency of any item in a stream using memory far below the number of distinct items.

Use \(d\) hash functions and a \(d\times w\) counter array. For each arrival, increment \(C[i, h_i(x)]\) for every \(i\). To query, report the minimum across rows:

\[\hat f(x)=\min_{i=1,\dots,d} C\big[i,\,h_i(x)\big].\]

Guarantee. With \(w=\lceil e/\varepsilon\rceil\) and \(d=\lceil\ln(1/\delta)\rceil\), the estimate satisfies \[f(x)\ \le\ \hat f(x)\ \le\ f(x)+\varepsilon N\] with probability at least \(1-\delta\), where \(N\) is the total stream length.

The estimate is never an underestimate, collisions only add counts, and taking the minimum over \(d\) independent hashes makes a large overestimate exponentially unlikely.

cms_new <- function(w, d, seed = 71) {
  set.seed(seed)
  list(C = matrix(0L, d, w), w = w, d = d,
       a = sample.int(2^30, d), b = sample.int(2^30, d))
}
cms_hash <- function(s, i, sk) {
  h <- sum(strtoi(charToRaw(as.character(s)), 16L) * seq_along(charToRaw(as.character(s))))
  (as.numeric(sk$a[i]) * h + sk$b[i]) %% 2147483647 %% sk$w + 1
}
cms_add <- function(sk, s) {
  for (i in seq_len(sk$d)) sk$C[i, cms_hash(s, i, sk)] <- sk$C[i, cms_hash(s, i, sk)] + 1L
  sk
}
cms_query <- function(sk, s)
  min(vapply(seq_len(sk$d), \(i) sk$C[i, cms_hash(s, i, sk)], numeric(1)))

set.seed(73)
N_cms <- 20000
# Zipf-distributed stream: a few very frequent items, a long tail
items <- sample(paste0("item", 1:3000), N_cms, TRUE,
                prob = 1 / (1:3000)^1.1)
truth <- table(items)

eps <- 0.001; delta <- 0.01
sk <- cms_new(w = ceiling(exp(1) / eps), d = ceiling(log(1 / delta)))
for (s in items) sk <- cms_add(sk, s)

probe <- names(sort(truth, decreasing = TRUE))[c(1, 2, 5, 20, 100, 500)]
data.frame(item = probe,
           true_count = as.integer(truth[probe]),
           sketch_estimate = vapply(probe, \(s) cms_query(sk, s), numeric(1)),
           error_bound = ceiling(eps * N_cms)) |>
  mutate(overestimate = sketch_estimate - true_count,
         within_bound = overestimate <= error_bound)

c(distinct_items = length(truth),
  exact_table_bytes = as.numeric(object.size(truth)),
  sketch_bytes = as.numeric(object.size(sk$C)),
  guarantee = sprintf("f(x) <= f_hat(x) <= f(x) + %d, w.p. >= %.2f",
                      ceiling(eps * N_cms), 1 - delta))

#>                                distinct_items 
#>                                        "2031" 
#>                             exact_table_bytes 
#>                                      "147616" 
#>                                  sketch_bytes 
#>                                       "54600" 
#>                                     guarantee 
#> "f(x) <= f_hat(x) <= f(x) + 20, w.p. >= 0.99"

Every estimate is at least the truth and within the bound, using a fixed-size array whose memory does not depend on the number of distinct items.

13.3 HyperLogLog

Problem. Count distinct elements in a stream. Exact counting requires \(\Theta(\text{distinct})\) memory; HyperLogLog uses a few kilobytes for cardinalities into the billions.

Hash each item to a uniform bit string. In a uniform stream, seeing a hash with \(\rho\) leading zeros suggests roughly \(2^{\rho}\) distinct items. Averaging \(\max\rho\) across \(m\) independent buckets, harmonically, to control variance — gives

\[\hat n=\alpha_m m^2\left(\sum_{j=1}^{m}2^{-M_j}\right)^{-1}, \qquad \text{relative error}\ \approx\ \frac{1.04}{\sqrt m}.\]

hll_estimate <- function(x, b = 12) {
  m <- 2^b
  h <- vapply(as.character(x), function(s) {
    r <- charToRaw(s); v <- sum(as.integer(r) * (31^(seq_along(r) - 1)))
    (v * 2654435761) %% 2^32
  }, numeric(1))
  idx <- (h %% m) + 1
  w <- floor(h / m)
  rho <- ifelse(w == 0, 32 - b, 32 - b - floor(log2(pmax(w, 1))))
  M <- rep(0, m)
  for (i in seq_along(idx)) M[idx[i]] <- max(M[idx[i]], rho[i])
  alpha <- if (m >= 128) 0.7213 / (1 + 1.079 / m) else 0.709
  est <- alpha * m^2 / sum(2^(-M))
  if (est <= 2.5 * m) { z <- sum(M == 0); if (z > 0) est <- m * log(m / z) }
  list(estimate = est, bytes = as.numeric(object.size(M)),
       theoretical_rel_error = 1.04 / sqrt(m))
}

set.seed(81)
hll_tab <- do.call(rbind, lapply(c(1e3, 5e3, 2e4, 5e4), function(nd) {
  x <- paste0("u", sample.int(nd * 3, nd * 4, TRUE))
  tru <- length(unique(x)); h <- hll_estimate(x, b = 12)
  data.frame(true_distinct = tru,
             hll_estimate = round(h$estimate),
             rel_error = round(abs(h$estimate - tru) / tru, 4),
             theoretical = round(h$theoretical_rel_error, 4),
             sketch_KB = round(h$bytes / 1024, 1),
             exact_KB = round(as.numeric(object.size(unique(x))) / 1024, 1))
}))
hll_tab

The sketch size is constant while the exact set grows linearly, and the observed relative errors sit near the theoretical \(1.04/\sqrt m\).

Sketch accuracy depends on both the memory budget and the stream length, which makes it a surface:

b_grid <- 6:14
n_grid <- c(1e3, 5e3, 1e4, 5e4, 1e5)

Zhll <- outer(b_grid, n_grid, Vectorize(function(b, nd) {
  set.seed(b + nd)
  x <- paste0("u", sample.int(nd * 3, min(nd * 2, 2e5), TRUE))
  tru <- length(unique(x))
  abs(hll_estimate(x, b = b)$estimate - tru) / tru
}))

plot_ly(x = n_grid, y = b_grid, z = Zhll, type = "surface",
        colorscale = "Inferno", reversescale = TRUE,
        colorbar = list(title = "Relative\nerror")) |>
  layout(title = "HyperLogLog relative error over memory budget and stream cardinality",
         scene = list(xaxis = list(title = "Distinct items", type = "log"),
                      yaxis = list(title = "log2(buckets m)"),
                      zaxis = list(title = "Relative error")))

Rotate along the bucket axis: error falls as \(1/\sqrt m\), doubling \(b\) quadruples \(m\) and halves the error. Rotate along the cardinality axis: the surface is nearly flat, which is the property that matters. The error does not grow with the stream.

14 Stream clustering

Clustering an unbounded stream cannot store the data, so the standard design is two-phase: an online phase maintains a compact summary in micro-clusters, and an offline phase reclusters those summaries on demand into macro-clusters.

library(stream)

set.seed(12345)
x_coords <- c(0.2, 0.3, 0.5, 0.8, 0.9)
y_coords <- c(0.8, 0.3, 0.7, 0.1, 0.5)
p_weight <- c(0.10, 0.30, 0.25, 0.20, 0.15)     # a proper probability vector

stream_5G <- DSD_Gaussians(k = 5, d = 2,
                           mu = cbind(x_coords, y_coords), p = p_weight)
stream_5G

#> Gaussian Mixture (d = 2, k = 5) 
#> Class: DSD_Gaussians, DSD_R, DSD

Note the weights: p is the mixture probability vector over the five components, so it must sum to 1. Supplying arbitrary positive numbers is silently renormalized by some generators and rejected by others, either way, the intended mixture is not what you get.

# D-Stream: grid-based density micro-clustering. A cell becomes a micro-cluster
# when its density exceeds Cm times the average cell density.
dstream <- DSC_DStream(gridsize = 0.1, Cm = 1.2)
update(dstream, stream_5G, n = 1000)

# Offline phase: recluster the micro-clusters into macro-clusters
km_macro <- DSC_Kmeans(k = 5)
recluster(km_macro, dstream)

c(micro_clusters = nclusters(dstream, type = "micro"),
  macro_clusters = nclusters(km_macro, type = "macro"),
  memory_KB = round(as.numeric(object.size(dstream)) / 1024, 1))

#> micro_clusters macro_clusters      memory_KB 
#>           17.0            5.0            3.2

plot(km_macro, stream_5G, type = "both",
     xlab = "X", ylab = "Y", main = "Micro-clusters (circles) and macro-clusters (crosses)")

Micro-clusters appear as circles sized by weight; macro-clusters as crosses. The summary occupies a few kilobytes regardless of how many points have streamed through, the defining property of the two-phase design.

14.1 Evaluating a stream clustering

Common misconception: “purity is a good measure of cluster quality.” Purity is \[\mathrm{Purity}=\frac1N\sum_{i=1}^{k}\max_j\big|c_i\cap t_j\big|,\] and it increases monotonically with \(k\): at \(k=N\) every cluster is a singleton and purity is exactly 1. It therefore cannot be used to choose the number of clusters, and it is not corrected for chance, a random partition into many clusters scores well.

The adjusted Rand index (Chapter 8, §8.4.2) has neither defect: it is 0 for a random partition and does not reward splitting.

set.seed(91)
n_p <- 600
truth_p <- rep(1:4, each = n_p / 4)
Xp <- cbind(rnorm(n_p, rep(c(0, 4, 0, 4), each = n_p / 4)),
            rnorm(n_p, rep(c(0, 0, 4, 4), each = n_p / 4)))

purity <- function(cl, tr)
  sum(vapply(unique(cl), \(i) max(table(tr[cl == i])), numeric(1))) / length(tr)

mono <- do.call(rbind, lapply(c(2, 4, 8, 16, 32, 64, 150), function(k) {
  cl <- kmeans(Xp, k, nstart = 10)$cluster
  data.frame(k = k, purity = purity(cl, truth_p),
             ARI = mclust::adjustedRandIndex(cl, truth_p))
}))
mono |> mutate(across(where(is.numeric), \(z) round(z, 4)))

mono |> pivot_longer(-k, names_to = "measure", values_to = "value") |>
  ggplot(aes(k, value, color = measure)) +
  geom_vline(xintercept = 4, linetype = "dashed", color = "grey45") +
  geom_line(linewidth = 1) + geom_point(size = 2.2) +
  scale_x_log10(breaks = mono$k) +
  scale_color_manual(values = c(purity = "#D8433B", ARI = "#3B7DD8")) +
  labs(title = "Purity rewards splitting; ARI does not",
       subtitle = "Dashed line: the true k = 4. Purity keeps climbing past it; ARI peaks there",
       x = "Number of clusters k (log scale)", y = NULL, color = NULL) +
  theme_dspa()

Purity climbs monotonically toward 1; ARI peaks at the true \(k=4\) and falls away. Only one of these can select a model.

14.2 Case study: knee-pain locations as a stream

The SOCR knee-pain dataset records \(x\)–\(y\) pain locations for over 8,000 patients, labelled by view (front/back × left/right).

library(rvest)
knee_raw <- dspa_try({
  rvest::read_html("https://wiki.socr.umich.edu/index.php/SOCR_Data_KneePainData_041409") |>
    html_nodes("table") |> _[[2]] |> html_table()
}, fallback = {
  set.seed(93)
  v <- rep(c("LF", "LB", "RF", "RB"), each = 500)
  data.frame(x = rnorm(2000, rep(c(30, 30, 70, 70), each = 500), 8),
             Y = rnorm(2000, rep(c(30, 70, 30, 70), each = 500), 8),
             View = v)
}, label = "SOCR knee-pain data")

rescale01 <- function(v) (v - min(v)) / (max(v) - min(v))

# data.frame(), not cbind(): a factor passed through cbind() is silently
# reduced to its integer codes and the level names are lost.
knee <- data.frame(x = rescale01(knee_raw$x),
                   y = rescale01(knee_raw$Y),
                   view = factor(knee_raw$View))
# 1. Rename the 'view' column to '.class' to support subsequent validation
names(knee)[names(knee) == "view"] <- ".class"

levels(knee$view); table(knee$view)

#> NULL

#> < table of extent 0 >

set.seed(1234)
knee <- knee[sample(nrow(knee)), ]     # views arrive in blocks; shuffle them

# The class label is passed via the `class=` argument, NOT as a third
# coordinate. Feeding it in as a feature and then evaluating against it would
# be circular (Chapter 8, Section 8.4.2).
# stream_knee <- DSD_Memory(knee[, c("x", "y")], class = knee$view, loop = TRUE)
# stream_knee <- DSD_Memory(knee[, c("x", "y", "view")], loop = TRUE)
stream_knee <- DSD_Memory(knee[, c("x", "y")], loop = TRUE)
stream_knee

#> Memorized Stream 
#> Class: DSD_Memory, DSD_R, DSD 
#> Contains 8666 data points - currently at position 1 - loop is TRUE

pts <- get_points(stream_knee, n = 5, info = TRUE)
pts

reset_stream(stream_knee, pos = 200)
head(get_points(stream_knee, n = 3, info = TRUE), 3)

reset_stream(stream_knee, pos = 1)
dsc_knee <- DSC_DStream(gridsize = 0.08, Cm = 0.6)
update(dsc_knee, stream_knee, n = 3000)

km_knee <- DSC_Kmeans(k = 4)             # four views, so four macro-clusters
recluster(km_knee, dsc_knee)

c(micro = nclusters(dsc_knee, type = "micro"),
  macro = nclusters(km_knee, type = "macro"))

#> micro macro 
#>    29     4

plot(km_knee, stream_knee, type = "both", xlim = c(0, 1), ylim = c(0, 1),
     main = "Knee-pain stream: micro-clusters and k-means macro-clusters")

DSPA3: Data Science and Predictive Analytics (UMich HS650)

Chapter 10: Data Formats, Streams, and Computational Performance

SOCR/MCAIM (Ivo Dinov)

August 28, 2026

How this chapter uses graphics