# A Simulation problem for a game

library(arm)

# (a)
pick=function()
{
	bag=c(rep('r',5),rep('g',5))
	return(sample(bag,5))
}
pick()
pick()
pick()

prize=function()
{
	balls=pick()
	nred=sum(balls=='r')
	if (nred==0 | nred==5) return(100-2)
	if (nred==1 | nred==4) return(5-2)
	return(-2)
}
prize()
prize()
prize()

# (b)
try1=replicate(10000,prize())
mean(try1)
sd(try1)

# (c)
try2=replicate(10000,prize())
count=0
for (i in 1:10000)
{
if(try2[i]>0) count=count+1
}
count/10000
mean(try2)

times=function(n)
{
	return(sum(replicate(n,prize())))
}
try3=replicate(10000,times(10))
count=0
for (i in 1:10000)
{
if(try3[i]>0) count=count+1
}
count/10000
# or
sum(try3>0)/10000
mean(try3)

# (d)
quantile(try3,seq(0.05,0.95,0.1))
range(try3)
histogram(try3,type='count',breaks=seq(range(try3)[1]-2.5,range(try3)[2]+2.5,5))

# (e)
X=replicate(10000,sum(pick()=='r')==4)
sum(X==T)/10000

# For reference
dhyper(0:5,5,5,5)
sum ( dhyper(0:5,5,5,5) * c( 100-2, 5-2, -2, -2, 5-2, 100-2 ) )



# Simulation for coefficients from a linear model

# (a)
chirps=read.table(file.choose(),header=T)
lm1=lm(chirp~temp, data=chirps)
display(lm1,digit=3)
lm0=lm(chirp~1, data=chirps)
anova(lm0,lm1)

# (b)
s=sim(lm1,1000)
print(apply(s$beta,2,quantile,probs=c(.025,.975)),digit=4)
0.216-2*0.039; 0.216+2*0.039

# (c)
set.seed(85)
sdsim=function(size) return(sd(sim(lm1,size)$beta[,2]))
sim10=(replicate(10,sdsim(10)))
sim100=(replicate(10,sdsim(100)))
sim1000=(replicate(10,sdsim(1000)))
sim10000=(replicate(10,sdsim(10000)))
mean(sim10); sd(sim10)
mean(sim100); sd(sim100)
mean(sim1000); sd(sim1000)
mean(sim10000); sd(sim10000)

# (d)
set.seed(85)
q0.01sim=function(size) return(quantile(sim(lm1,size)$beta[,2],0.01))
sim10=(replicate(10,q0.01sim(10)))
sim100=(replicate(10,q0.01sim(100)))
sim1000=(replicate(10,q0.01sim(1000)))
sim10000=(replicate(10,q0.01sim(10000)))
mean(sim10); sd(sim10)
mean(sim100); sd(sim100)
mean(sim1000); sd(sim1000)
mean(sim10000); sd(sim10000)

# (e)
pre=s$beta%*%as.matrix(c(1,70))+rnorm(1000,0,s$sigma)
print(quantile(pre,c(.025,.5,.975)),digit=3)
print(predict(lm1,data.frame(temp=70)),digit=3)
display(lm1)
print(14.5-2*.98,digit=3); print(14.5+2*.98,digit=3)