# lets set n to some value. Play around with various sizes of n

n <- 30

#
# Lets generate samples from the standard normal distribution and
# examine normal probability plots
#
#


par(mfrow=c(3,3))
for (i in 1:9){
  x <- rnorm(n)
  qqnorm(x)
  qqline(x)
}



#
# Less to do with lab 1 but lets do normal probability plots
# of non normal data

# exponential data

par(mfrow=c(3,3))
for (i in 1:9){
 x <- rexp(n)
 qqnorm(x)
 qqline(x)
}


# t distribution with df=9

par(mfrow=c(3,3))
for (i in 1:9){
 x <- rt(n,9)
 qqnorm(x)
 qqline(x)
}



# cauchy distribution 

par(mfrow=c(3,3))
for (i in 1:9){
 x <- rcauchy(n,9)
 qqnorm(x)
 qqline(x)
}


# gamma distribution shape=2, rate=3

par(mfrow=c(3,3))
for (i in 1:9){
 x <- rgamma(n,2,3)
 qqnorm(x)
 qqline(x)
}


##
## we can construct ourselves a qq function for
## checking if data is gamma with certain
## shape and rate parameters
##
##

gamma.qqplot <- function(x,shape=2,rate=3){

  theoretical.probs <- seq(1:length(x))/(length(x)+1)
  theoretical.quantiles <- qgamma(theoretical.probs,shape,rate)

  plot(theoretical.quantiles, sort(x),xlab="Theoretical Quantiles",ylab="Sample Quantiles",main="Gamma QQ-plot")
  # add a line through the 25% and 75% quantiles (in same manner as qqline)
 
  y <- quantile(x,prob=c(0.25,0.75))
  x <- qgamma(c(0.25,0.75),shape,rate)
  slope <- diff(y)/diff(x)
  int <- y[1] - slope*x[1]
  abline(int,slope)
}


##
## Now what happens with normals
##

par(mfrow=c(3,3))
for (i in 1:9){
  x <- rnorm(n)
  gamma.qqplot(x)
}


##
## How about with gammas?
##
##

par(mfrow=c(3,3))
for (i in 1:9){
 x <- rgamma(n,2,3)
 gamma.qqplot(x)
}