library(arm)
# (a)
pick=function()
{
bag=c(rep('r',5),rep('g',5))
return(sample(bag,5))
}
bag=c(rep('r',5),rep('g',5))
bag
pick()
pick()
pick()
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)
}
balls=pick()
balls
balls=='r'
sum(balls=='r')
nred=sum(balls=='r')
nred==0
nred==0 | nred==5
nred==2
prize()
prize()
prize()
prize()
prize()
prize()
prize()
try1=replicate(10000,prize())
mean(try1)
sd(try1)
sum(try1)
quantile(try1)
try2=replicate(10000,prize())
count=0
for (i in 1:10000)
{
if(try2[i]>0) count=count+1
}
count/10000
count
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
time(10)
time(10)
times(10)
times(10)
times(10)
times(10)
times(10)
times(10)
times(10)
times(10)
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))
sum(replicate(10000,sum(pick()=='r')==4)==T)/10000
sum(pick()=='r')
sum(pick()=='r')==4
replicate(10000,sum(pick()=='r')==4)
sum(replicate(10000,sum(pick()=='r')==4)==T)
sum(replicate(10000,sum(pick()=='r')==4)==T)/10000
chirps=read.table(file.choose(),header=T)
chirps
lm1=lm(chirp~temp, data=chirps)
display(lm1,digit=3)
lm0=lm(chirp~1, data=chirps)
anova(lm0,lm1)
s=sim(lm1,1000)
s
s$beta
s$sigma
print(apply(s$beta,2,quantile,probs=c(.025,.975)),digit=4)
0.216-2*0.039; 0.216+2*0.039
apply(s$beta,2,quantile,probs=c(.025,.975))
display(lm1)
# (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)
sdsim=function(size) return(sd(sim(lm1,size)$beta[,2]))
sdsim=function(size) { return(sd(sim(lm1,size)$beta[,2])) }
# (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)
predict(lm1,data.frame(temp=70))
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)
pre=s$beta%*%as.matrix(c(1,70))+rnorm(1000,0,s$sigma)
print(quantile(pre,c(.025,.5,.975)),digit=3)
pre