##############################################################
#  read.table example                                        #
##############################################################

treatment.df <- read.table(file="treatment_data.txt")
# You may need to specify the exact path to the file, e.g.
# treatment.df <- read.table(file="C:/Documents and Settings/Administrator/Desktop/Rcourse/treatment_data.txt")

# Note also that blank line and lines starting with # are ignored by default
treatment.df <- read.table(file="treatment_data.txt", header=T)


group
attach(treatment.df)
group

edit(treatment.df)

s.rate = n.resp/n.group
cbind(n.group,n.resp,s.rate)

plot(duration, s.rate)
plot(duration, s.rate,type="n")
points(duration[type=="A"], s.rate[type=="A"], pch="A")
points(duration[type=="B"], s.rate[type=="B"], pch="B", col="red")

detach(treatment.df)





##############################################################
# Simulation of some distributions                           #
# and Kolmogorov-Smirnof gof test                            #
##############################################################


### Exponential
lambda = 2.0
y1 = rexp(200, rate = lambda)
par(mfrow=c(1,2))
hist(y1, col="cyan",main="Histogram of Y1 ~ Exp(2)")
boxplot(y1, horizontal=T, col="cyan",main="Boxplot of Y1")
descriptives <- list(summary(y1), var(y1));
descriptives


### Gamma
alpha = 3.0; beta = 2.0
y2 = rgamma(200, shape = alpha, rate = beta)
hist(y2, col="cyan",main="Histogram of Y2 ~ Gamma(3,2)")
boxplot(y2, horizontal=T, col="cyan",main="Boxplot of Y2")
descriptives <- list(summary(y2), var(y2));
descriptives


### LN
simulate.ln.f <- function(n,mu,sigma2){
  y3 = exp(rnorm(n,mean=mu, sd=sqrt(sigma2)))
  #
  par(mfrow=c(1,2))
  hist(y3, col="cyan",main=paste("Histogram of Y3 ~ LN(", mu, ",", sigma2,")"))
  boxplot(y3, horizontal=T, col="cyan",main="Boxplot of Y3")
  #
  descriptives <- list(summary(y3), var(y3));
  #
  return(descriptives)
}

simulate.ln.f(n=200, mu=0, sigma2=0.1)


### Weibull
##  Function simulating from Weibull distn
weib.r <- function(n, nu, lambda){
	u = runif(n,0,1)
	r = ( - log(1-u)/lambda )^(1/nu)
	return(r)
}


simulate.weib.f <- function(n, nu, lambda){
  y4 = weib.r(n,nu,lambda)
  #
  par(mfrow=c(1,2))
  hist(y4, col="cyan",main=paste("Histogram of Y4 ~ Weib(", nu, ",", lambda,")"))
  boxplot(y4, horizontal=T, col="cyan",main="Boxplot of Y4")
  #
  descriptives <- list(summary(y4), var(y4));
  #
  return(list(y4,descriptives))
}

simulate.weib.f(n=200, nu=2, lambda=0.5)



###
### Perform KS test and plot empirical v theoretical CDF
###


#  Function returning cdf for Weibull distn
weib.cdf <- function(q, nu, lambda){
	cdf =  1- exp(-lambda*q^nu)
	return(cdf)
}


ks.weib.f <- function(data,nu,lambda){
  #  Perform test
  ks <- ks.test(data,weib.cdf,nu,lambda)
  #  Plot ecdf and cdf
  grid.x = seq(min(data), max(data), length=100)
  par(mfrow=c(1,1))
  plot(grid.x,weib.cdf(grid.x,nu,lambda),type="l",col="red", ylim=c(0,1))
  s = c(1:length(data))
  lines(sort(data),s/length(data), type="s")
  title(main="Empirical v theoretical CDF")
  legend("bottomright", legend=c("cdf","ecdf"),col=c("red","black"),lty=c(1,1))
  #
  return(ks)
}


weib.data = simulate.weib.f(n=200, nu=2, lambda=0.5)[[1]]
#ks.weib.f(weib.data, nu=2, lambda=0.5)
ks.weib.f(weib.data, nu=2, lambda=0.4)