##########################################################
# Author: Oliver Elison Timm
# Date: 2014-01-27
# 
# Purpose: Demonstration of Probability laws
# Frequentist's approach to measure probability 
#
##########################################################


# use the random generator to draw numbers from 1 to 6
# ( our 6-sided rolling die)
# but now we repeat this experiment over different
# trial sizes and see how the randomness
# and sample size affect our estimated probabilities
nside<-6

ntrial<-c(10,20,30,40,50,100,200,500,1000,10000)
iexp<-length(ntrial)
# iexp is the number of ntrial numbers 
# for each of the ntrial-repeated experiments we will store
# the estimated probability for the die events 1-6
# here we use a 2-d matrix with columns for the
# events 1,2,3,4,5,6
# and rows for the experiements with different ntrial
Pest<-matrix(NA,iexp,nside)
for (k in 1:iexp) {
  count<-rep(0,nside)
  n<-ntrial[k]
  for (i in 1:n) {
    # we only pick the first element of this random sequence 
    res<-sample(nside)
    event<-res[1]
    count[event]<-count[event]+1
  }
  print("----------------------------------")
  print(paste("estimated with ntrials=",n))
  print("frequency of events:")
  print(seq(1,nside,1))
  print(count/n)
  Pest[k,1:nside]<-count/n
}
print(paste("theoretical probability P(e): 1/6",1/6))


# A simple plot to see how close our estimates are
plot(c(1,6),c(0.0,0.35),typ='n',xlab='event',ylab='est. prob.')

# a short version for seq(1,nside,1)
x<-1:nside
for (k in 1:iexp) {
  # pch selects a different symbol
  points(x,Pest[k,],pch=k)
}
# adding the theoretical values as line
lines(x,rep(1/6,nside),col=2)
# adding a legend
x<-1.15
y<-0.36
label<-as.character(ntrial)
lchar<-1:iexp
legend(x,y,label,pch=lchar)
dev.copy2pdf(file="figures/prob_rolldie_trial_stat.pdf")