# Today we explore the following event:

# Event: A randomly selected group of 29 people contain at least one  pair of people that share a birthday in common. 

#Below (Offerings) we see what we would offer (as a market) for the bet of receiving 1 dollar if  the above Event occurs.

#  The first column of "Offerings"  is the price you agree to pay someone to take this bet, while  the second column is the price you agree to sell this bet.  The third column is the second minus the first.  If this were negative you would be offering free money in your market!   As a market maker, this is a bad idea. The size of the difference  measures our confidence in the bet's true "Fair Price".    If the true Fair price is in your interval, then in the long run (if people keep placing bets) you will make money; while if the Fair price is outside your buy and sell range, then eventually someone will be sure to make money off you!  This bet's fair price is the probability of this event occurring and we will explore this idea further.

Guesses<-c(.3,.4,1/23,2/23,.1,.2,.06,.14,0.03,0.07,.07,.1,.22,.42,.5,.65,.25,.75,.28,.41,.033,.055,.01,.08,21/365,23/365,.25,.35,.08,0.12,(1.16)*10^(-59),(1.18)*10^(-59),0.05,0.15)
dim(Guesses)<-c(2,length(Guesses)/2)
G<-rbind(Guesses,Guesses[2,]-Guesses[1,])
O<-order(G[3,])
Offerings<-round(aperm(G[,O],c(2,1)),3)
Offerings

#Notice many people are incompatible! What is the true fair price?  Below we see what happened to the average number of successes in many experiments.  Namely, we simulated the experiment assuming that a "random" person's Birthday  is equal likely to be one of the 365 possibilities (we simplified away leap years). 

set.seed(c(12072004,8023333271),kind ="Marsaglia-Multicarry")
ClassSize=23
Trials=1000
BDs<-ceiling(365*runif(ClassSize*Trials))
dim(BDs)<-c(ClassSize,Trials)
sortBDs<-apply(BDs,2,sort)
Test<-sortBDs[2:ClassSize,]-sortBDs[1:(ClassSize-1),]
Ex<-apply(Test,2,min)
CA=cumsum(Ex==0)/(1:Trials)
plot(1:Trials,CA,type='l',main="Long Range Tendency of Birthday Bet",xlab="Number of Trials",ylab="Fraction of Occurrences")
sum(Ex==0)/Trials

# We found in class that after 10^3 experiments the rather unexpected 0.521  (unexpected since only one confident person had the upper end of their range greater than 0.50).  After 10^5 we found 0.50641 On my linux box (using the built in generator and seed) I ran it for 10^6  experiments and found  0.507211.     How close are we?

# Well in class we found the true answer is:

Days=365
Size=23
1-(1/365^365)*(factorial(Days)/factorial(Days-Size))

#Oooops!  What happened! (Overflow).  We better be a little smarter. 

X<-(Days-Size+1):Days
Y<-rep(Days,Size)
Prob=1-prod(X/Y)
Prob

#Hence the true answer is very nearly:

0.5072972

# A graph might be nice as well

Days=365
Sizes=1:100
Probs<-c()
for (Size in Sizes){
X<-(Days-Size+1):Days
Y<-rep(Days,Size)
Probs=c(Probs,1-prod(X/Y))
}

plot(Sizes,Probs,type="l",main="Birthday coincidence probabilities", ylab="Birthday coincidence probabilities",xlab="Sample Size")


#################################
Class Birthdays
#################################

#Here was our actual class.  Did we share any birthdays? (Yes!) How many? (Two pairs! How likely is that? Try to compute and if you cant then simulate!)

Class<-list(c("Jan",23),c("Sep",9),c("Oct",19),c("Oct",16),c("Oct",9),c("Oct",26),c("Apr",4),c("May",18),c("Apr",16),c("Jun",17),c("Dec",15),c("Sep",12),c("Mar",27),c("Dec",1),c("May",22),c("Dec",15),c("Apr",24),c("Mar",27),c("May",12),c("Sep",14),c("May",29))

C<-Dates(Class)
SC<-sort(C)
Class[C==349]


######## Dates is a function that converts a list in this form into our funny number birthdays.

makenum<-function(L){
 {if (L=="Jan")
 x=1}
  {if (L=="Feb")
 x=2}
  {if (L=="Feb")
 x=2} 
 {if (L=="Mar")
 x=3}
  {if (L=="Apr")
 x=4} 
 {if (L=="May")
 x=5} 
 {if (L=="Jun")
 x=6}
  {if (L=="Jul")
 x=7}
  {if (L=="Aug")
 x=8} 
 {if (L=="Sep")
 x=9}
 {if (L=="Oct")
 x=10}
 {if (L=="Nov")
 x=11}
 {if (L=="Dec")
 x=12}
 x
}

Con<-function(m,d){
ML<-c(31,28,31,30,31,30,31,31,30,31,30,31)
t=-ML[1]
for (i in 1:m){
t=ML[i]+t
}
t<-t+d
t
}

Date<-function(L){
Mon=L[1]
day=L[2]
Con(makenum(Mon),as.numeric(day))
}

Date(c("Jan",10))

Dates<-function(L){
N=length(L)
H<-c()
for (i in 1:N){
H<-c(H,Date(L[[i]]))
}
H
}



l<-length(Class)
C<-Dates(Class)
SC<-sort(C)
T<-SC[2:l]-SC[1:(l-1)]
sum(T==0)
Class[C==349]