Search code examples
rprobability

Illustrating the LLN (Law of Large Numbers)


I have to illustrate the Law of Large Numbers through simulations in R. More precisely. I would like to illustrate that the cumulative distribution function of the mean,

enter image description here

converges to the function f given by f(x) = 0 if x ≤ μ and f(x) = 1 if x > μ. In my case, I have to use a dice. That is, each Xi is the uniformly distributed on {1,2,3,4,5,6}, so μ = 3.5.
Using R, I have tried to proceed in the following way:

n <- 100
N <- 10000
mu <- 3.5
for(j in 1:N)
{
V[j] <- sum(sample(1:6), n, replace = TRUE);
}
f <- function(x)
{
if (x<=3.5)
{
y <-0
}
else
{
y <- 1
}
y
}
Vf <- Vectorize(f, "x")

So my idea was to compare the cumulative distribution function of the mean with the function f using a plot. How can I implement it in R properly. So I have to plot the cumulative distribution function and the function f in one plot.


Solution

  • You can simulate dice-rolls like this

    set.seed(1)
    n.rolls <- 100
    dicerolls <- sample(1:6, n.rolls, replace=TRUE)
    
    mean(dicerolls)
    

    As for the rest of your question I'm afraid I'd need some further explanation. Maybe you can draw an image of what kind of plot you want?

    If this is homework you should tag your question accordingly, and read the info for the tag.
    As you can see this site doesn't support MathJax/LaTeX equation mark-up. If you want to include equations you can do it through something like codecogs.

    Maybe it's something like this you're thinking of?

    dicerolls <- function(rolls=2, reps=10^4) {
        mean.per.replicate <- replicate(reps, mean(sample(1:6, rolls, replace=TRUE)))
    }
    
    set.seed(1)
    dice.seq <- c(1:6, 20, 100)
    
    opar <- par(no.readonly=TRUE)
    par(mar=c(2, 2.5, 1, 0.1), mfrow=c(length(dice.seq), 2), 
      cex=0.5, mgp=c(1.5, 0.5, 0))
    
    for (i in dice.seq) {
    
        hist(dicerolls(i), breaks=50, col="darkgrey", 
          xlim=c(1, 6), ylim=c(0, 3), freq=FALSE, main="", xlab="")
        legend("topleft", paste(i, "dice"), bty="n")
    
        plot(ecdf(dicerolls(i)), xlim=c(1, 6), main="", frame.plot=FALSE)
    }
    
    par(opar)
    

    enter image description here