#=====================================================
# Simulation of the law of large numbers.
#=====================================================

#-----------------------------------------------------
# Function to generate sequence of successes and
# failures.  Arguments:  number of trials is n, success
# probability on any given trial is p.  Function returns
# vector x of 0's and 1's.
#-----------------------------------------------------
outcomes=function(n,p) {
  u=runif(n)
  x=1*(u<=p)   # R converts logical vector to numeric 0's
               #  and 1's in attempting to multiply by 1.
  return(x)
}

#-----------------------------------------------------
# Simulate every number of trials from 5 to max.n.
# Store proportions of successes in p.obs.
#-----------------------------------------------------
p=.26          # Probability of success.
max.n=500      # Maximum number of trials.
n=5:max.n
p.obs=numeric(length(n))
for (i in 1:length(n)) {
  p.obs[i]=sum(outcomes(n[i],p))/n[i]
}
#-----------------------------------------------------
# Plot p.obs versus n.
#-----------------------------------------------------
plot(n,p.obs,type="l",lty=1,xlab="number of trials",
  ylab="proportion of successes")
abline(h=p,lty=2)