SMaC: Statistics, Math, and Computing

APSTA-GE 2006: Applied Statistics for Social Science Research

Eric.Novik@nyu.edu | Summer 2026 | Session 3

Session 3 Outline

  • Transforming data for plotting
  • Deriving the free-fall equation with antiderivatives
  • Some rules of integration
  • Evaluating integrals numerically
  • The waiting time distribution

\[ \DeclareMathOperator{\E}{\mathbb{E}} \DeclareMathOperator{\P}{\mathbb{P}} \DeclareMathOperator{\V}{\mathbb{V}} \DeclareMathOperator{\L}{\mathscr{L}} \DeclareMathOperator{\I}{\text{I}} \]

Data transformations

  • We saw how annual and continuous compounding compare at a fixed interest rate
  • How does the asset value vary with the interest rate?
  • We investigate with a common simulate-pivot-plot pattern

Map Function

“The purrr::map* functions transform their input by applying a function to each element of a list or atomic vector and returning an object of the same length as the input.”

  • Generate 10 vectors of 100 uniform random values. The first vector has min = 1, the second has min = 2, and so on through min = 10; all have max = 15
  • Then compute the mean of each of the 10 vectors
library(purrr)

set.seed(2006)
x <- 1:10
samples <- x |>
  map(\(minimum) runif(n = 100, min = minimum, max = 15))

means <- samples |>
  map_dbl(mean)
  • Run the code and inspect samples and means. How do the two objects differ, and what did you expect to see?

Pivot Functions

pivot_longer() “lengthens” data, increasing the number of rows and decreasing the number of columns. The inverse transformation is pivot_wider().

library(tidyr)
head(iris)
  Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1          5.1         3.5          1.4         0.2  setosa
2          4.9         3.0          1.4         0.2  setosa
3          4.7         3.2          1.3         0.2  setosa
4          4.6         3.1          1.5         0.2  setosa
5          5.0         3.6          1.4         0.2  setosa
6          5.4         3.9          1.7         0.4  setosa
iris_long <- iris |> pivot_longer(!Species, names_to = "length_width", values_to = "measure")
head(iris_long)
# A tibble: 6 × 3
  Species length_width measure
  <fct>   <chr>          <dbl>
1 setosa  Sepal.Length     5.1
2 setosa  Sepal.Width      3.5
3 setosa  Petal.Length     1.4
4 setosa  Petal.Width      0.2
5 setosa  Sepal.Length     4.9
6 setosa  Sepal.Width      3

Data Transformations

library(purrr)
library(tidyr)
library(dplyr)
library(ggplot2)

rates <- seq(0.05, 0.20, length = 10)
P <- 100
time <- seq(1, 50, length = 50)

Pe <- function(A, r, t) A * exp(r * t)
d <- time |>
  map(\(x) Pe(A = P, r = rates, t = x))
class(d)
[1] "list"
d[[1]][1:4]
[1] 105.1271 106.8939 108.6904 110.5171
names(d) <- as.character(time)
d <- as_tibble(d)
1 2 3
105.1271 110.5171 116.1834
106.8939 114.2631 122.1403
108.6904 118.1360 128.4025
110.5171 122.1403 134.9859
112.3745 126.2802 141.9068
114.2631 130.5605 149.1825
116.1834 134.9859 156.8312
118.1360 139.5612 164.8721
120.1215 144.2917 173.3253
122.1403 149.1825 182.2119
  • Modern R usage offers many shortcuts, but those may be confusing to beginners.
  • In particular, R loops have mostly been replaced with map() functions.
  • We recommend purrr::map() functions instead of R’s *apply().
rates <- seq(0.05, 0.20, length = 10)
P <- 100
time <- seq(1, 50, length = 50)

Pe <- function(A, r, t) A * exp(r * t)
time |> map(\(x) Pe(A = P, r = rates, t = x))

# above is a shortcut for
map(time, function(x) Pe(A = P, r = rates, t = x))

# and the above is a shortcut for the following loop
l <- list()
for (i in seq_along(time)) {
  l[[i]] <- Pe(A = P, r = rates, t = time[i])
}
library(tidyr)
# add rates as a column
d <- d %>% mutate(rate =
            round(rates, 2) |>
              as.character())

# convert from wide format to long
d <- d %>%
  pivot_longer(!rate,
               names_to = "year",
               values_to = "value")

d$year <- as.numeric(d$year)
  • Here is a cheat sheet explaining tidyr functions.
  • And here is a much shorter version of the exercise.
rate year value
0.05 1 105.1271
0.05 2 110.5171
0.05 3 116.1834
0.05 4 122.1403
0.05 5 128.4025
0.05 6 134.9859
0.05 7 141.9068
0.05 8 149.1825
0.05 9 156.8312
0.05 10 164.8721 
p <- ggplot(d, aes(year, value))
p + geom_line(aes(color = rate), linewidth = 0.2) +
  scale_y_continuous(labels =
          scales::dollar_format()) +
  xlab("Time (years)") +
  ylab("Asset value") +
  ggtitle("Growth of $100 at different interest rates")

d |>
  filter(year == max(year)) |>
  select(rate, y50 = value) |>
  knitr::kable()
rate y50
0.05 1218.249
0.07 2803.162
0.08 6450.009
0.1 14841.316
0.12 34149.510
0.13 78577.199
0.15 180804.241
0.17 416026.201
0.18 957266.257
0.2 2202646.579

RStudio cheat sheets

Your Turn

  • Run install.packages("HistData")
  • Then run library(HistData)
  • Take a look at the Arbuthnot dataset: ?Arbuthnot
Year Males Females Plague Mortality Ratio Total
1629 5218 4683 0 8771 1.114243 9.901
1630 4858 4457 1317 10554 1.089971 9.315
1631 4422 4102 274 8562 1.078011 8.524
1632 4994 4590 8 9535 1.088017 9.584
1633 5158 4839 0 8393 1.065923 9.997
1634 5035 4820 1 10400 1.044606 9.855
  • Use these tools to produce a plot that looks something like this
  • Bonus: compute and plot the proportion of christenings that were female, \(\text{Females}/(\text{Males}+\text{Females})\)

Integral Calculus

  • Integration describes accumulated change and plays a central role in statistics
  • For continuous distributions, integrals compute probabilities and expectations
  • Across Bayesian and frequentist statistics, integration is used to normalize densities and marginalize over unknown quantities
  • Derivatives are often used for optimization, including maximum-likelihood and maximum-a-posteriori estimation

  • Unknown quantities are often called parameters, such as the local gravitational acceleration \(g\) in our Moon-drop model

Intuition Behind Integration

  • A definite integral is a continuous analog of summation: it measures accumulated signed change over an interval
  • Geometrically, it is the signed area between a one-dimensional function \(f\) and the horizontal axis
  • An indefinite integral represents a family of antiderivatives: \(\int f(x)\,dx = F(x)+C\) when \(F'(x)=f(x)\)
  • The Fundamental Theorem of Calculus connects the two: \(\int_a^b f(x)\,dx = F(b)-F(a)\)
  • Many definite integrals are evaluated numerically, but it helps to understand what the computer is approximating

Image source: Calculus Volume 1

Some Common Integrals

OpenStax: Here is a more complete list

Techniques of Integration

  • Definite integration, like differentiation, is linear:

\[ \begin{aligned} \int_a^b [f(x)+g(x)]\,dx &= \int_a^b f(x)\,dx + \int_a^b g(x)\,dx, \\ \int_a^b c f(x)\,dx &= c\int_a^b f(x)\,dx. \end{aligned} \]

  • Unlike differentiation, integration has no universal procedure that produces a closed-form antiderivative for every function

  • Many common integrals have closed-form analytical solutions, but many important integrals in statistics do not

  • For well-behaved integrands in one or two dimensions, numerical quadrature is often practical

  • Higher-dimensional problems may use Monte Carlo, importance sampling, quasi-Monte Carlo, or Markov chain Monte Carlo (MCMC)

  • For simple integrals, we can sometimes find a closed-form solution using substitution or integration by parts

Playing with Integrals

  • Given our intuition for integrals as signed areas, let’s see how to compute them analytically and numerically
  • The numerical techniques shown here work best in low dimensions; high-dimensional integration requires other methods
  • Suppose we want to evaluate the integral \(\int_{1}^{3} x \sin(x^2)\,dx\)

Constructing a Riemann Sum

riemann <- function(f, lower, upper, n = 100) {
  edges <- seq(lower, upper, length.out = n + 1)
  total <- 0

  for (i in seq_len(n)) {
    width <- edges[i + 1] - edges[i]
    height <- f(edges[i + 1]) # right endpoint
    total <- total + width * height
  }

  total
}

Notice that the function takes another function as an argument. Functions that do this are called higher-order functions.

Source: OpenStax Calculus Volume 1

Evaluating the Integral

f <- function(x) x * sin(x^2)
x <- seq(0, pi, length.out = 100)

riemann(f, 1, 3, n = 200)
[1] 0.7275414
# compute using R's integrate function
integrate(f, 1, 3)
0.7257163 with absolute error < 1.3e-09

Integrating Analytically

  • Many integrals cannot be evaluated analytically, but we can find an antiderivative of \(x\sin(x^2)\)
  • Let \(u=x^2\). Then \(du=2x\,dx\), so \(x\,dx=\tfrac{1}{2}du\)

\[ \begin{aligned} \int x\sin(x^2)\,dx &= \tfrac{1}{2}\int \sin(u)\,du \\ &= -\tfrac{1}{2}\cos(u)+C \\ &= -\tfrac{1}{2}\cos(x^2)+C. \end{aligned} \]

Comparing the Results

  • Using R’s integrate function:
integrate(f, 1, 3)
0.7257163 with absolute error < 1.3e-09
  • Using the analytical solution
f1 <- function(x) -1/2 * cos(x^2)

f1(3) - f1(1)
[1] 0.7257163
  • The universe is in balance!

Analytical Integration on the Computer

  • When in doubt, you can always try WolframAlpha

  • The Python library SymPy can be used through the R package caracas

  • A computer algebra system returns one antiderivative and usually omits the arbitrary constant \(+C\)

library(caracas); library(stringr)
add_align <- function(latex) {
  str_c("\\begin{align} ", latex, " \\end{align}")
}
add_int <- function(latex) {
  str_c("\\int ", latex, "\\, dx")
}
x <- symbol('x'); f <- x^2 / sqrt(x^2 + 4)
tex(f) %>% add_int() %>% str_c(" =") %>% add_align() %>% cat()

\[\begin{align} \int \frac{x^{2}}{\sqrt{x^{2} + 4}}\, dx = \end{align}\]

int(f, x) %>% tex() %>% add_align() %>% cat()

\[\begin{align} \frac{x \sqrt{x^{2} + 4}}{2} - 2 \operatorname{asinh}{\left(\frac{x}{2} \right)} \end{align}\]

Your Turn: Waiting Time

  • The exponential distribution is commonly used to model waiting times between events that occur at a constant average rate \(\lambda>0\)
  • Its probability density function (PDF) is

\[ f(x)= \begin{cases} \lambda e^{-\lambda x}, & x\ge 0, \\ 0, & x<0. \end{cases} \]

  • Every PDF must integrate to 1. Verify that

\[ \int_0^\infty \lambda e^{-\lambda x}\,dx = 1. \]

  • The exponential distribution is memoryless: conditional on having already waited \(s\), the probability of waiting at least another \(t\) does not depend on \(s\)

\[ \P(X>s+t\mid X>s)=\P(X>t)=e^{-\lambda t}. \]

Integration by Parts

The product rule gives a useful integration identity:

\[ \begin{aligned} (fg)' &= f'g + fg', \\ \int f(x)g'(x)\,dx &= f(x)g(x)-\int f'(x)g(x)\,dx, \\ \int u\,dv &= uv-\int v\,du. \end{aligned} \]

For example, integration by parts gives the mean waiting time:

\[ \begin{aligned} \E[X] &= \int_0^\infty x\lambda e^{-\lambda x}\,dx \\ &= \left[-xe^{-\lambda x}\right]_0^\infty + \int_0^\infty e^{-\lambda x}\,dx \\ &= \frac{1}{\lambda}. \end{aligned} \]

Exponential Growth (again)

The exponential-growth function used earlier solves the following differential equation. Assume \(y(t)>0\):

\[ \frac{\text{d}[y(t)]}{\text{d}t} = k \cdot y(t) \]

We can now solve it:

\[ \begin{align*} \frac{1}{y} \, \text{d}y &= k \, \text{d}t \\ \int \frac{1}{y} \, \text{d}y &= \int k \, \text{d}t \\ \log(y) &= k \cdot t + C \\ y(t) &= y_0 \cdot e^{kt}, \, y_0 = e^C \end{align*} \]

Deriving the Free-Fall Equation

  • Let \(x(t)\) be the downward distance fallen from the release point, so downward is positive
  • Treat the gravitational acceleration \(g\) as constant near the surface
  • The ball has constant mass \(m\), and the only force acting on it is gravity, \(F=mg\)

Newton’s second law relates force to the rate of change of momentum:

\[ \begin{aligned} F &= ma \\ &= \frac{d}{dt}\left(m\frac{d[x(t)]}{dt}\right) = \frac{d}{dt}\left(mv(t)\right)\\ &= mg. \end{aligned} \]

Because \(m\) is constant, take it outside the derivative and cancel it from both sides:

\[ m\frac{d^2[x(t)]}{dt^2}=mg \quad\Longrightarrow\quad \frac{d^2[x(t)]}{dt^2}=g. \]

Deriving the Free-Fall Equation

Integrate twice with respect to \(t\):

\[ \begin{aligned} \int \frac{d^2[x(t)]}{dt^2}\,dt &= \int g\,dt, \\ \frac{d[x(t)]}{dt} &= gt+C_1, \\[0.6em] \int \frac{d[x(t)]}{dt}\,dt &= \int (gt+C_1)\,dt, \\ x(t) &= \tfrac{1}{2}gt^2+C_1t+C_2. \end{aligned} \]

The ball is released from rest, and distance is measured from the release point, so

\[ \left.\frac{d[x(t)]}{dt}\right|_{t=0}=0 \Longrightarrow C_1=0, \qquad x(0)=0 \Longrightarrow C_2=0. \]

Therefore,

\[ \boxed{x(t)=\tfrac{1}{2}gt^2}. \]

Free Fall on Earth’s Moon

Using the accepted value \(g_{\text{Moon}}=1.625\ \text{m/s}^2\), the ball reaches the bottom when \(x(t)=100\) m:

\[ \begin{aligned} 100 &= \tfrac{1}{2}g_{\text{Moon}}t^2, \\ t_{\text{impact}} &= \sqrt{\frac{2(100)}{g_{\text{Moon}}}} \\ &\approx 11.1\ \text{s}. \end{aligned} \]

The free-fall model describes the motion only until impact.

Homework

  • Take a look at the iris dataset (?iris)
  • Compute the overall means of Sepal.Length, Sepal.Width, Petal.Length, and Petal.Width
  • Compute these means by Species (hint: use group_by() and summarise() from dplyr)
  • Produce a plot that looks like this:

  • Produce a plot that looks like this (see geom_density() and facet_wrap()):

  • Compute the following integral and show the steps:

\[ \int 2x \cos(x^2)\, dx \]

  • Evaluate this integral on paper from \(0\) to \(2\pi\), and use R’s integrate() function to validate your answer.