Simulation - part 2 Stat 133 Gaston Sanchez Department of - PowerPoint PPT Presentation

Simulation - part 2 Stat 133 Gaston Sanchez Department of Statistics, UC–Berkeley gastonsanchez.com github.com/gastonstat Course web: gastonsanchez.com/stat133

Random Numbers in R 2

Generating Random Numbers Generation of random numbers is at the heart of many statistical methods 3

Use of Random Numbers Some uses of random numbers ◮ Sampling procedures ◮ Simulation studies of stochastic processes ◮ Analytically intractable mathematical expressions ◮ Simulation of a population distribution by resampling from a given sample from that population ◮ More general: Simulation, Monte Carlo, Resampling 4

Random Samples ◮ Many statistical methods rely on random samples: – Sampling techniques – Design of experiments – Surveys ◮ Hence, we need a source of random numbers ◮ Before computers, statisticians used tables of random numbers ◮ Now we use computers to generate “random” numbers ◮ The random sampling required in most analyses is usually done by the computer 5

Generating Random Numbers ◮ We cannot generate truly random numbers on a computer ◮ Instead, we generate pseudo-random numbers ◮ i.e. numbers that have the appearance of random numbers ◮ they seem to be randomly drawn from some known distribution ◮ There are many methods that have been proposed to generate pseudo-random numbers 6

Generating Random Numbers A very important advantage of using pseudo-random numbers is that, because they are deterministic, they can be reproduced (i.e. repeated) 7

Multiple Recursion ◮ Generate a sequence of numbers in a manner that appears to be random ◮ Use a deterministic generator that yields numbers recursively (in a fixed sequence) ◮ The previous k numbers determine the next one x i = f ( x i − 1 , . . . , x i − k ) 8

Simple Congruential Generator ◮ Congruential generators were the first reasonable class of pseudo-random number generators ◮ The congruential method uses modular arithmetic to generate “random” numbers 9

Ingredients ◮ An integer m ◮ An integer a such that a < m ◮ A starting integer x 0 (a.k.a. seed ) 10

Simple Congruential Generator The first number is obtained as: x 1 = ( a × x 0 ) mod m The rest of the pseudorandom numbers are generated as: x n +1 = ( a × x n ) mod m 11

Simple Congruential Generator For example if we take m = 64 , and a = 3 , then for x 0 = 17 we have: x 1 = (3 × 17) mod 64 = 51 x 2 = (3 × 51) mod 64 = 25 x 3 = (3 × 25) mod 64 = 11 x 4 = (3 × 11) mod 64 = 33 x 5 = (3 × 33) mod 64 = 35 x 6 = (3 × 35) mod 64 = 41 . . . 12

Congruential algorithm a <- 3; m <- 64; seed <- 17 x <- numeric(20) x[1] <- (a * seed) %% m for (i in 2:20) { x[i] <- (a * x[i-1]) %% m } x[1:16]; x[17:20] ## [1] 51 25 11 33 35 41 59 49 19 57 43 1 3 9 27 17 ## [1] 51 25 11 33 13

Congruential algorithm cong <- function(n, a = 69069, m = 2^32, seed = 17) { x <- numeric(n) x[1] <- (a * seed) %% m for (i in 2:n) { x[i] <- (a * x[i-1]) %% m } x } y <- cong(20, a = 3, m = 64, seed = 17) 14

Congruential algorithm plot(y[1:(n-1)], y[2:n], pch = 19, xlab = 'current value', ylab = 'next value') 60 ● ● 50 ● ● ● 40 ● next value ● ● ● 30 ● ● ● 20 ● ● 10 ● ● ● ● ● 0 0 10 20 30 40 50 60 current value 15

cong(n, a = 69069, m = 2^32, seed = 17) 4e+09 3e+09 next value 2e+09 1e+09 0e+00 0e+00 1e+09 2e+09 3e+09 4e+09 current value 16

Ingredients ◮ m is the modulus ( 0 < m ) ◮ a is the multiplier ( 0 < a < m ) ◮ x 0 is the seed ( 0 ≤ x 0 ≤ m ) The period length of a Random Generator Number is at most m , and for some choices of a much less than that. 17

About the seed ◮ You can reproduce your simulation results by controlling the seed ◮ set.seed() allows to do this ◮ Every time you perform simulation studies indicate the random number generator and the seed that was used 18

About the seed # set the seed set.seed(69069) runif(3) # call the uniform RNG ## [1] 0.1648855 0.9564664 0.3345479 runif(3) # call runif() again ## [1] 0.01109596 0.18654873 0.94657805 # set the seed back to 69069 set.seed(69069) runif(3) ## [1] 0.1648855 0.9564664 0.3345479 19

Simple Congruential Generator We can use a congruential algorithm to generate uniform numbers We’ll describe one of the simplest methods for simulating independent uniform random variables on the interaval [0, 1] 20

Generating Uniform Numbers The generator proceeds as follows x 1 = ( ax 0 ) mod m u 1 = x 1 /m u 1 is the first pseudo-random number, taking some value between 0 and 1. 21

Ingredients The second pseudorandom number is obtained as: x 2 = ( ax 1 ) mod m u 2 = x 2 /m u 2 is another pseudorandom number. 22

Generating Uniform Numbers ◮ If m and a are chosen properly, it is difficult to predict the value u 2 given u 1 ◮ For most practical purposes u 2 is approximately independent of u 1 23

Simple Congruential Generator For example if we take m = 7 , and a = 3 , then for x 0 = 2 we have: x 1 = (3 × 2) mod 7 = 6 , u 1 = 0 . 857 x 2 = (3 × 6) mod 7 = 4 , u 2 = 0 . 571 x 3 = (3 × 4) mod 7 = 5 , u 3 = 0 . 714 x 4 = (3 × 5) mod 7 = 1 , u 4 = 0 . 143 x 5 = (3 × 1) mod 7 = 3 , u 5 = 0 . 429 x 6 = (3 × 3) mod 7 = 2 , u 6 = 0 . 286 . . . 24

Random Numbers in R R uses a pseudo random number generator ◮ It starts with a seed and an algorithm (i.e. function) ◮ The seed is plugged into the algorithm and a number is returned ◮ That number is then plugged into the algorithm and the next number is created ◮ The algorithm is such that the produced numbers behave like random numbers 25

RNG functions in R Function Description Uniform runif() rbinom() Binomial Multinomial rmultinom() Negative binomial rnbinom() Poisson rpois() Normal rnorm() rbeta() Beta Gamma rgamma() rchisq() Chi-squared Cauchy rcauchy() See more info: ?Distributions 26

Random Number Functions runif(n, min = 0, max = 1) sample from the uniform distribution on the interval (0,1) The chance the value drawn is: ◮ between 0 and 1/3 has chance 1/3 ◮ between 1/3 and 1/2 has chance 1/6 ◮ between 9/10 and 1 has chance 1/10 27

Random Number Functions rnorm(n, mean = 0, sd = 1) sample from the normal distribution with center = mean and spread = sd 28

Random Number Functions rnorm(n, mean = 0, sd = 1) sample from the normal distribution with center = mean and spread = sd rbinom(n, size, prob) sample from the binomial distribution with number of trials = size and chance of success = prob 28

More Simulations 29

Simulation Probability examples For instance ◮ Simulate flipping a coin ◮ Simulate rolling a die ◮ Simulate drawing a card from a deck ◮ Simulate a probability experiment with balls in an urn ◮ Simulate the “Monty Hall Problem” 30

Flipping a Coin 31

Simulating a Coin How to simulate tossing a coin? 32

Simulating a coin One way to simulate a coin coin <- c("heads", "tails") 33

Simulating a coin One way to simulate a coin coin <- c("heads", "tails") One way to simulate flipping a coin sample(coin, size = 1) ## [1] "heads" 33

Probability Probability allows us to quantify statements about the chance of an event taking place For example, flip a fair coin ◮ What’s the chance it lands heads? ◮ Flip it 4 times, what proportion of heads do you expect? ◮ Will you get exactly that proportion? ◮ What happens when you flip the coin 1000 times? 34

Simulating a Coin When you flip a coin ◮ it may land heads ◮ it may land tails ◮ with what probability it lands heads? ◮ If it is a fair coin: p = 0 . 5 35

Simulating a Coin Tossing a coin can be modeled with a random variable following a Bernoulli distribution: ◮ heads ( X = 1 ) with probability p ◮ tails ( X = 0 ) with probability q = 1 − p The Bernoulli distribution is a special case of the Binomial distribution: B (1 , p ) 36

Simulating tossing a coin Tossing a coin simulated with a binomial distribution: # binomial distribution generator rbinom(n = 1, size = 1, prob = 0.5) ## [1] 0 37

Flipping a Coin function Function that simulates flipping a coin # flipping coin function coin <- function(prob = 0.5) { rbinom(n = 1, size = 1, prob = prob) } coin() ## [1] 1 38

Flipping a Coin function It’s better if we assign labels # flipping coin function coin <- function(prob = 0.5) { out <- rbinom(n = 1, size = 1, prob = prob) ifelse(out, "heads", "tails") } coin() ## [1] "tails" 39

Flipping a Coin function # 10 flips for (i in 1:10) { print(coin()) } ## [1] "heads" ## [1] "tails" ## [1] "tails" ## [1] "heads" ## [1] "tails" ## [1] "tails" ## [1] "tails" ## [1] "tails" ## [1] "tails" ## [1] "heads" 40

4 Flips In 4 flips ◮ Possible outputs: – HHHH, THHH, HTHH, HHTH, HHHT, ... ◮ we can get 0, 1, 2, 3, 4 heads ◮ so the proportion of heads can be: 0, 0.25, 0.5, 0.75, 1 ◮ we expect the proportion to be 0.5 ◮ but a proportion of 0.25 is also possible 41