Outline

◼
  • Uniform RNGs
    • ◼
  • What is a pseudo-random number
    • ◼
  • Plugging one into the framework
    • ◼
  • Continuous Distributions
    • ◼
  • Functional Transformations
    • ◼
  • Inversion
    • ◼
  • Acceptance Rejection
    • ◼
  • Discrete Distributions

    • Randomness

    ◼
  • What it means for something to be random is a bit philosophical but a way to look at it is a complicated process where we don’t have access to the details and we have averaged them out.
    • ◼
  • Computers are (we wish) deterministic so generating a (pseudo) random variable with a computer involves transformations that garble the inputs enough to lose track of where we started.
    • ◼
  • Samplers from a uniform (pseudo) random distribution are the fundamental building blocks we use to build samplers from more complicated distributions, we want them too look as random as possible balancing other factors such as speed, memory, and simplicity.
    • Histogram[RandomReal[{0,1},10^6]]

    Uniform RNGs in the WL

    ◼
  • RandomInteger
    — random integer or array of integers
    • ◼
  • RandomReal
    — random real or array of reals
    • ◼
  • RandomComplex
    — random complex numbers
    • ◼
  • RandomChoice
    — random choices from a list
    • ◼
  • RandomSample
    — random samples or permutation of a list
    • ◼
  • RandomPrime
    — random prime number
    • ◼
  • RandomColor
    ▫
    RandomImage
    ▫
    RandomGraph
    ▫
    RandomPermutation
    ▫
    RandomEntity
    ▫
    RandomWord

    • Seeding and Precision

    ◼
  • Two important aspects of the framework are:
    • ◼
  • We can specify the seeding to get reproducible results
    • ◼
  • We can get arbitrary precision numbers with more random digits than those allowed by hardware
    • SeedRandom[2];RandomReal[10]​​SeedRandom[2];RandomReal[10]
    7.2224
    7.2224
    RandomReal[1]​​Precision[%]​​RandomReal[1,WorkingPrecision->24]​​Precision[%]
    0.196193
    MachinePrecision
    0.711858589835955524055212
    24.

    Uniform RNGs - Plug your own

    ◼
  • The RNG framework allows for one to plug-in a custom uniform RNG
    • ◼
  • https://reference.wolfram.com/language/tutorial/RandomNumberGeneration.html#Defining_Your_Own_Generator
    • ◼
  • Let’s have a look at PCG which is an RNG with very strong statistical properties among other nice things:
    • ◼
  • This is code adapted from: https://mathematica.stackexchange.com/a/91731/6581
    • PCG[bw_,mul_,inc_,const_,state_]["GenerateBits"[bits_]]:=​​Module[{bot,ns,res,rot,xsh,maskBW,mask2BW}​​,​​mask2BW=BitShiftLeft[1,2bw]-1;​​maskBW=BitShiftLeft[1,bw]-1;​​{xsh,bot,rot}=const;​​​​ns=BitAnd[mulstate+inc,mask2BW];​​​​res=BitShiftRight[BitXor[state,BitShiftRight[state,xsh]],bot];​​res=BitRotateRight[BitAnd[res,maskBW],BitShiftRight[state,rot],bw];​​​​{res,PCG[bw,mul,inc,const,ns]}​​]
    BlockRandom[​​SeedRandom[1,Method->PCG];​​data=RandomReal[1,10^4]​​];​​Histogram[data]

    Non-uniform RNGs

    ◼
  • To sample from non-uniform RNGs (distributions) we can use RandomVariate
    • ◼
  • RandomVariate works for both continuous and discrete distributions
    • data=RandomVariate[NormalDistribution[],10^6];
    data=RandomVariate[PoissonDistribution[3],10^6];
    Poisson PDF

    Functional Transformation

    ◼
  • The next step is to generate samples from a non-uniform distribution using a uniform RNG.
    • ◼
  • Random variables can be fully described by their CDF
    • F
    X
    (x)=P(X≤x)
    ◼
  • This formula allows us to construct random variables by applying a function to other random variables
    • ◼
  • We can use TransformedDistribution to figure out the distribution of these transformed random variables

    • Inversion Formula

    ◼
  • What we want is to generate samples from X   having access to samples of U  
    • ◼
  • What is g?
    • ◼
  • For a uniform random variable U the CDF takes the form:

    • Inversion Formula

    ◼
  • Let’s construct samples from the normal distribution using its inverse CDF

    • Box-Muller Transform

    ◼
  • As we saw in the last slide evaluating certain inverse CDFs can be very expensive.
    • ◼
  • There are all sort of clever transformations to sample from important distributions.
    • ◼
  • Let’s look at one which takes 2 uniform random number and mixes them to produce 2 Normally distributed random numbers.
    • ◼
  • This function can be inverted using elementary special functions, which often have hardware support.
    • ◼
  • We can then use Pythagoras Theorem

    • Acceptance-Rejection

    ◼
  • Another very general method uses 2 uniform random numbers per sample of X:
    • ◼
  • one to propose a possible sample from a “hat” distribution H.
    • ◼
  • Let’s sample from a Trapezoidal Distribution with a uniform hat:

    • Acceptance-Rejection Algorithm

    ProbabilityDistribution and Acceptance-Rejection

    ◼
  • ProbabilityDistribution we saw earlier is a general container for user defined distributions.
    • ◼
  • It has several sampling methods available, including acceptance-rejection
    • ◼
  • We normalised the density earlier but it is not necessary for sampling and expensive

    • Discrete Distributions-Poisson

    ◼
  • Let’s talk a little a bit discrete distributions
    • ◼
  • Here’s an algorithm by Knuth that samples from a Poisson distribution by halting when a threshold is reached.
    • ◼
  • In our system we use a more complicated algorithm, which switches between inversion and an Acceptance Rejection algorithm by Ahrens and Dieter, which uses discrete normal for a hat function
    • ◼
  • The hat is above the Poisson density on one (known) side of the range and below for the other

    • Acceptance-Rejection Closing Remarks

    ◼
  • I showed a toy implementation of the algorithm, but different facets of Acceptance-Rejection are used in most sampling algorithms
    • ◼
  • Acceptance-Rejection algorithms can get very clever and be very efficient.
    • ◼
  • It becomes expensive when the hat gets “far” away from the density.
    • ◼
  • This happens inevitably at higher dimensions: “Curse of dimensionality”

    • Closing Remarks

    ◼
  • Huge subject - Just scratches the surface
    • ◼
  • Coming up talks:
    • ◼
  • Part 2 - RandomProcesses
    • ◼
  • Part 3 - Point Processes
    • ◼
  • Part 4 - Random Matrices?
    • ◼
  • Deliberately avoided talking about MCMC
    • ◼
  • Higher dimensions: Point Processes

    • CITE THIS NOTEBOOK

    Wolfram R&D LIVE: Probability and Statistics - Random Sampling​
    by Eduardo Serna​
    Wolfram Community, STAFF PICKS, June 27, 2023
    ​https://community.wolfram.com/groups/-/m/t/2946101