a little adjustment. # LLN epsilon = 0.1 m = 100 # repetition to get the probability n = 1000 percent = 1:n index = 1:n # the number of r.s. to get xbar for (i in 1:n){ Count = 0 for (j in 1:m){ z = rnorm(i) y = mean(z) if (abs(y) > epsilon) Count = Count + 1 } percent[i] = Count/m } plot(index, percent, type='p', main = 'Approximate Probability of exceeding epsilon') percent[n]