setwd("C:/Documents and Settings/William Christensen/My Documents/stat230")

9 + 5
x <- c(5, 10, 20, 40, 800)
x
mean(x)   # this is a "comment"
median(x)  # the median is a better measure of spread 
           # when extreme skewness is exhibited
sd(x) 
help(sd)
help(c)

normfunction <- function(x,mu=0,sigma=1)  # creating a function
{
  1/sqrt(2*pi*sigma^2) * exp(-1/(2*sigma^2) * (x-mu)^2)
}

normfunction(0)   # something like this already exists as a pre-existing function
                  #  in R: try dnorm(0)


xpts <- seq(-10,10,.01)
xpts
plot(xpts,normfunction(xpts,-5,1),type="l",col="blue",lwd=2, xlab="x", ylab="f(x)")
lines(xpts,normfunction(xpts,0,4),col="red",lwd=2)
lines(xpts,normfunction(xpts,5,3),col="darkgreen",lwd=2)
legend(4,.35,c("N(-5, 1)","N(0, 4)","N(5, 3)"),col=c("blue","red","darkgreen"),lwd=2)

plot(xpts,dnorm(xpts))

xpts <- seq(-3.5,3.5,.01)
plot(xpts,normfunction(xpts,0,1),type="l",lwd=2, xlab="x", ylab="f(x)")
segments(x0=seq(-3,3,1),y0=0,y1=normfunction(seq(-3,3,1)))
abline(h=0)


### 
### One-sample t-test
###
x <- rnorm(20,0.3,1)   ## Generate 20 random numbers from the N(0.3,1) distribution
                       ##    where 0.3 is the mean and 1 is the STD. DEV. (not variance)
mean(x)
sd(x)
n <- length(x)
hist(x)

## Test Ho: mu = 0 vs. Ha: mu > 0 (greater than 0)
t <- (mean(x) - 0) / (sd(x)/sqrt(n))
t
pval <- 1-pt(t,n-1)
pval

## Test Ho: mu = 0 vs. Ha: mu != 0 (not equal to 0)
pval <- 2*(1-pt(abs(t),n-1))
pval

# or 
help(t.test)
# EASY WAY OF DOING 1-SAMPLE T-TEST
t.test(x, alternative="greater")

## 98% confidence interval for mu
lowercl <- mean(x) - qt(.99,n-1) * sd(x) / sqrt(n)
uppercl <- mean(x) + qt(.99,n-1) * sd(x) / sqrt(n)
c(lowercl,uppercl)
# or...EASY WAY:
t.test(x, alternative="two.sided", conf.level=.98)


###
### Working problems 41 and 42 from practice exam
###
spl2 <- ((35-1)*5.05^2 + (56-1)*4.13^2) / (35+56-2)
spl2
sqrt(spl2) 

(6.21-3.66) - qt(.95,89) * sqrt( (1/35 + 1/56) * spl2 )
(6.21-3.66) + qt(.95,89) * sqrt( (1/35 + 1/56) * spl2 )

t <- (6.21-3.66) / sqrt( (1/35 + 1/56) * spl2 )
t

mydf <- 35+56-2

2*(1-pt(abs(t),mydf))




### 
### Two-sample t-test: Shoe Material Example
###

## TYPE DATA IN...  ##
A <- c(13.2, 8.2, 10.9, 14.3, 10.7, 6.6, 9.5, 10.8, 8.8, 13.3)
B <- c(14.0,8.8,11.2,14.2,11.8,6.4,9.8,11.3,9.3,13.6)

## ...OR PASTE IN FROM FILE ##
# Read Method #1 (preferred for reproducibility): pasting data in:
shoes <- read.table(header=TRUE, text="
Material Wear
A 13.2
A 8.2
A 10.9
A 14.3
A 10.7
A 6.6
A 9.5
A 10.8
A 8.8
A 13.3
B 6.4
B 9.8
B 11.3
B 9.3
B 13.6
B 14
B 8.8
B 11.2
B 14.2
B 11.8
")

# OR Read Method #2: reading from a file stored locally
shoes <- read.table("shoes.txt", header=TRUE)

# OR Read Method #3: reading from a webpage (DON'T USE THIS ONE)
# The next 3 lines are needed if you want to read files directly from https (i.e., encrypted) webpages
install.packages("RCurl")   # Only run this once on a given machine...downloads the add-on package called "RCurl"
library("RCurl")   # this loads the add-on package called "RCurl" 
options(download.file.method="curl")  #Alters the way encrypted webpages are read, enabling the ability to read directly
shoes <- read.table("https://tofu.byu.edu/stat230/shoes.txt", header=TRUE)

A <- shoes[1:10,2]  ## rows 1 through 10, column 2
B <- shoes[11:20,2]

boxplot(cbind(A,B))  ## "column bind"
par(mfrow=c(1,2))
hist(cbind(A),nclass=8)
hist(cbind(B),nclass=8)

mean(A)
mean(B)
sd(A)
sd(B)
spl2 <- ((10 - 1)*var(A) + (10 - 1)*var(B)) / (10 + 10 - 2)
df <- 10 + 10 - 2

t <- (mean(A) - mean(B))/ sqrt( (1/10 + 1/10) * spl2 )
t
pval <- 2*(1 - pt(abs(t), df))
pval


# EASY WAY OF DOING 2-SAMPLE T-TEST
t.test(A,B, alternative="two.sided", var.equal=TRUE)




###
### Power curve for shoe data (two-sample test)
###

deltas <- seq(-6,0,.1)    #sequence of numbers from -6 to 0 by 0.1
powers <- rep(NA,length(deltas))    # create a vector that repeats NA for a 
                                    # number of times equal to the length of deltas
criticalvalue <- qt(.05, 18)* sqrt( (1/10 + 1/10) * spl2 ) 
     ## Finds the 5th percentile of the t distribution with 18 degrees of freedom and
     ##   multiplies it by the standard error of xbar1 - xbar2 to find out what 
     ##   value of xbar1 - xbar2 is required to have a t statistic in the lowest
     ##   5% of the distribution
criticalvalue

for (i in 1:length(deltas))
{
  tvalue <- (criticalvalue - deltas[i]) / sqrt( (1/10 + 1/10) * spl2 )
     ##  Finding how many standard errors away from (mu1 - mu2) is the criticalvalue
     ##     for different values of delta = mu1 - mu2.  
  powers[i] <- pt( tvalue , 18)
     ##  This is the probability of getting a value less than tvalue when drawing 
     ##     from a t distribution with 18 degrees of freedom
}

plot(deltas, powers, type="l", ylim=c(0,1), main="Power curve")
abline(h=0.05, col="red")
text( -4, .07, "0.05", col="red")

#plot(deltas, 1-powers, type="l", ylim=c(0,1), main="OC curve")






##  Paired comparison with Shoe Wear Data
shoespaired <- read.table(header=TRUE, text="
Material  Wear
A	13.2
A	8.2
A	10.9
A	14.3
A	10.7
A	6.6
A	9.5
A	10.8
A	8.8
A	13.3
B	14
B	8.8
B	11.2
B	14.2
B	11.8
B	6.4
B	9.8
B	11.3
B	9.3
B	13.6                    
")

# HARD WAY OF DOING PAIRED T TEST #
A <- shoespaired[1:10,2]  ## rows 1 through 10, column 2
B <- shoespaired[11:20,2]
d <- A - B
mean(d)
sd(d)
n <- length(A)

t <- mean(d) / (sd(d) / sqrt(n))
t

pval <- 2 * pt(-3.349, 10-1 )
pval

## OR ##
# EASY WAY OF DOING PAIRED T TEST #
t.test(A, B, paired=TRUE, alternative="two.sided")