#############################################
# First things
#############################################

motif()                  # starts a graphics window
par(mar=c(4,4,1,1))      # margins for plot (not necessary)
options(digits=3)        # sets precision to 3 digits (default is 7)
help.start()             # starts up an interactive help window
                          
#############################################
# Quick examples
#############################################
x _ geyser$waiting      
y _ geyser$duration
plot(x,y,xlab="WAITING",ylab="DURATION",pch="+",col=2)

for (i in 1:10){
  cat(i,"! = ",gamma(i+1),"\n")
}

#############################################
# Extra problem #2 on HW2
#############################################
# to find the cutoffs for the binomial & neg binomial:                     
pbinom(0:12,12,0.5)     # look at "binomial" in the help window
pnbinom(0:20,9,0.5)
rbind(0:12, pbinom(0:12,12,0.5)) # to save you the counting

#############################################
# Problem 3.3
#############################################
# To simulate from the difference of two t-dist r.v.'s:

n <- 1000
mu1 <- 0; s1 <- 1; mu2 <- 2; s2 <- 1; nu1 <- 18; nu2 <- 18
   # dummy values for a quick example -- see the problem for
   # the correct values
x1 <- rchisq(n,nu1)
x2 <- rchisq(n,nu2)
z1 <- rnorm(n)
z2 <- rnorm(n)
mut <- mu1+s1*z1*sqrt(nu1)/sqrt(x1)
muc <- mu2+s21*z2*sqrt(nu2)/sqrt(x2)
mut.minus.muc <- mut-muc
hist(mut.minus.muc,xlab="D = MU-t - MU-c",prob=T,ylab="P(D|y)")

# See the appendix, p481, for the simulation of the t. See the Splus
# help function for the functions rchisq(.) & rnorm(.)