p.empirical.distribution <- function(samplesize){

  simulations <- 1000
  true.p <- 0.01

  ## simulate observation of fever
  fever <- rbinom(simulations, samplesize, true.p)

  ## estimate p from the observation
  p.hat <- fever/samplesize

  ## histogram is saved as an object... 
  histogram <-  hist(p.hat,
                     breaks = seq(0,0.2,0.0005),
                     plot = FALSE
                     )
  
  ## ..and then frequency is changed to percent
  histogram$counts <- histogram$counts/simulations*100 

  ## now we can plot the histogram...
  plot(histogram,
       main = "",
       xlab = "Estimate",
       ylab = "Percent")    

  ## ...and add a box around the plot, just to make it look nice.
  box()
}

## draw empirical distribution of p for 4 differenct samplesizes
par(mfrow=c(2,2))
p.empirical.distribution(10)
p.empirical.distribution(100)
p.empirical.distribution(1000)
p.empirical.distribution(10000)