R makes it easy to generate random numbers from a wide variety of distributions with a consistent interface. For example, runif, rnorm, and rpois generate random numbers from uniform, normal, and Poisson distributions, respectively.

> x = runif(n = 1000, min = 4.6, max = 9.3)
> min(x)
[1] 4.60374
> max(x)
[1] 9.288063
> 
> y = rnorm(n = 1000, mean = 81, sd = 9)
> mean(y)
[1] 81.2548
> sd(y)
[1] 9.179427
> 
> z = rpois(n = 1000, lambda = 7.11)
> mean(z)
[1] 7.187
> var(z)
[1] 6.954986

The Racket math library provides much of the same functionality as the base R functions, but I first want to work through a couple of examples without using the math library because I think it is useful for better understanding Racket.

Let's create a Racket function that is similar to runif in R.

#lang racket/base

(define (runif n [min 0] [max 1])
  (for/list ([i (in-range n)])
    (+ (* (- max min) (random)) min)))

The random function draws a random number from a uniform distribution between 0 and 1. A list of n random samples is drawn by looping through a sequence with for/list. To draw a random number within a specified interval, we multiply the random number by the specified range (- max min) and shift by the min.

> (define unif (runif 1000 -7 33))
> (apply min unif)  ; check that min is about -7
-6.987705647350004
> (apply max unif)  ; check that max is about 33
32.97250528872971

The need to use apply to calculate the min of a list was a surprise to me. The min function takes any number of real numbers as arguments but not a list. apply uses the contents of the list as the arguments for min.

> (min 3 1 2)
1
> (min '(3 1 2))
min: contract violation
  expected: real?
  given: '(3 1 2)

Another wrinkle is introduced if you want to generate a vector (with for/vector) of random numbers rather than a list.

(define (runif-vec n [min 0] [max 1])
  (for/vector ([i (in-range n)])
    (+ (* (- max min) (random)) min)))

We can't use apply with a vector.

> (define unif-vec (runif-vec 1000 0.5 1.5))
> (apply min unif-vec)
apply: contract violation
  expected: list?
  given: '#(1.3465423772765357 1.0623864994795045 1.3437265664094888 1.1048441535810904 0.7861090566755002 0.6470981723620602 0.5546597487687198 0.552946824350613 0.5251621516034304 1.3568612624022975 1.1217164560959263 1.3334770245391927 0.890117942389224...
  argument position: 2nd
  other arguments...:

StackOverflow provides a solution based on for/fold.

> (for/fold ([m +inf.0]) ([x (in-vector unif-vec)]) (min m x))
0.5014616191163698
> (for/fold ([m -inf.0]) ([x (in-vector unif-vec)]) (max m x))
1.4994848021056595

In for/fold, m is an accumulator that is initialized as positive (+inf.0) or negative infinity (-inf.0) for min and max, respectively. When iterating through unif-vec, m is compared to x and the value of m is updated with the output of min or max.

Random values can also be drawn from a normal distribution using draws from the standard uniform distribution with random and a formula from Rosetta Code

> (define (rnorm n [mean 0] [sd 1])
   (for/list ([i (in-range n)])
    (+ mean (* (sqrt (* -2 (log (random)))) (cos (* 2 pi (random))) sd))))
> (define norm (rnorm 1000 2 0.5))
> (require math)      ; need math library for mean and standard deviation
> (mean norm)         ; check that mean is about 2
1.9797112583083123
> (stddev norm)       ; check that sd is about 0.5
0.5022756727111091

Now, let's simplify our life and use the functions provided by the math library.

(require math) 

(define (rnorm-2 n [mean 0] [sd 1])
  (if (and (real? mean) (real? sd))
      (flnormal-sample (real->double-flonum mean) (real->double-flonum sd) n)
      (error "mean and sd arguments must be real numbers")))

The function, flnormal-sample, requires that the mean and standard deviation are passed as floating-point numbers (and returns flvector). In rnorm-2, we check if the mean and sd arguments are real numbers and, if so, convert them to floating-point numbers (real->double-flonum) before calling flnormal-sample, which is otherwise similar to R's runif function.

> (define norm-2 (rnorm-2 1000 2 0.5))
> (mean norm-2)
2.024086073688603
> (stddev norm-2)
0.5039448247564015

The math library contains similar sample functions for numerous types of distributions, including uniform, Poisson, beta, binomial, gamma, and more.

Finally, let's take a look at the statistics object provided by the math library.

> (define s (update-statistics* empty-statistics norm-2))
> (statistics-mean s)
1.9947322187693184
> (statistics-stddev s)
0.4974799838849746

empty-statistics is an empty statistics object. update-statistics* populates the statistics object by iterating through norm-2 to compute summary statistics, which are extracted with accessor functions, e.g., statistics-min, statistics-mean, statistics-skewness.

The question is why would you use the more verbose statistics object. One reason is that statistics-min and statistics-max are actually less verbose than using for/fold with min and max. Also, if you are computing several summary statistics on a large sequence, then update-statistics* is faster than the other assorted functions.

> (define norm-3 (rnorm-2 10000000 1000 100))
> (time
   (mean norm-3)
   (stddev norm-3)
   (variance norm-3)
   (for/fold ([m +inf.0]) ([x (in-flvector norm-3)]) (min m x))
   (for/fold ([m -inf.0]) ([x (in-flvector norm-3)]) (max m x)))
cpu time: 72618 real time: 108130 gc time: 34828
> (time
   (define s2 (update-statistics* empty-statistics norm-3))
   (statistics-mean s2)
   (statistics-stddev s2)
   (statistics-variance s2)
   (statistics-min s2)
   (statistics-max s2))
cpu time: 12478 real time: 22668 gc time: 5512