#Problem 50 Chapter 9

#manganese<-read.table("U:/Math 152/MANGNESE.txt",header=T)
#manganese<-read.table("C:/Documents and Settings/Mike/Desktop/Math 152/MANGNESE.txt",header=T)
manganese<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/MANGNESE.txt",header=T)
names(manganese)
rm(d)
rm(c1)
rm(c2)
rm(c3)
rm(c4)
rm(c5)

attach(manganese)


 plot(d,c1)

 plot(d,c2)

 plot(d,c3)

 plot(d,c4)

 plot(d,c5)

 plot(d,(c1+c2+c3+c4+c5)/5)

 

mdata<-c(c1,c2,c3,c4,c5)
 hist(mdata)
  length(mdata)
#[1] 120


plot(qnorm((1/121)*(1:120)),mdata[order(mdata)])
cor(qnorm((1/121)*(1:120)),mdata[order(mdata)])
#[1] 0.9953098

f<-function(x){l<-rnorm(120)
cor(qnorm((1/121)*(1:120)),l[order(l)])
}


d<-sapply(1:10000,f)
d[order(d)[50]]
d[order(d)[100]]
d[order(d)[250]]
d[order(d)[500]]
d[order(d)[1000]]
d[order(d)[2000]]
d[order(d)[2500]]
d[order(d)[5000]]
d[order(d)[7500]]
d[order(d)[9000]]
d[order(d)[9500]]
d[order(d)[9750]]
d[order(d)[9900]]
d[order(d)[9950]]

 shapiro.test(mdata)

 #       Shapiro-Wilk normality test

#data:  mdata 
#W = 0.9892, p-value = 0.4658

 sqrt(.9892)
#[1] 0.9945853
 
#Use the following url  http://www.jstor.org/view/00359254/
# then select volume 31, 1982, no.2 pp. 115-124  
# to see a paper by Royston on "extending the Wilk shapiro test to large samples".
#You must be on a campus computer to do this. (or otherwise have a jstor account)

mean(mdata)
#[1] 1.409667
 sd(mdata)
#[1] 0.1522047

 hist(mdata,prob=T)
 curve(dnorm(x,m=1.409667,sd=.1522047),add=T)

phat<- pnorm((1 + .1*(1:8) -mean(mdata))/sd(mdata)) - pnorm((1 + .1*(0:7) -mean(mdata))/sd(mdata))
nhat<-length(mdata)*phat

n<-function(j){length(mdata[1 + (j-1)*.1 <mdata& mdata<=1+ j*.1])}
#Note that it is important to use "<=" in the above line since there are several data points which lie on the boundary of 
#the bins.
n<-sapply(1:8,n)

barplot(sqrt(n)-sqrt(nhat))
barplot((n-nhat)/sqrt(nhat))

# The fluctuations are all less than two standard deviations from what would be expected for normality.
# But the rootogram shows that the deviations may not be at random, there being too many observations in both
#the upper and lower tails. Given our preliminary plotting, we might suspect that the cause is the  evident drift downward in time of the observed
#values. The downward drift must also contribute to the variance of the observations.

#----------------------------------------------------------------------
#Problem 63 chapter 9

#whales<-read.table("U:/Math 152/whales.txt",header=T)
whales<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/whales.txt",header=T)

attach(whales)

names(whales)

#first with the mle estimates

alphahatmle<- 1.595414

lambdahatmle<-2.632701

length(times)
#[1] 210

hist(times,prob=T)
 curve(dgamma(x,shape=alphahatmle,rate=lambdahatmle),add=T)


plot(qgamma((1/211)*(1:210),shape=alphahatmle,rate=lambdahatmle),times[order(times)])

id<-function(x){x}

curve(id,add=T)

#now with the method of moments estimates

lambdahatm<- 1.318771

alphahatm<- 0.799174

hist(times,prob=T)
curve(dgamma(x,shape=alphahatm,rate=lambdahatm),add=T)


plot(qgamma((1/211)*(1:210),shape=alphahatm,rate=lambdahatm),times[order(times)])
curve(id,add=T)

#recall how it looked for the rain fall data

illrain<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/ILLRAIN.txt",header=F,na.strings="*")


attach(illrain)

rain<-c(V1,V2,V3,V4,V5)

rain<-rain[!is.na(rain)]

rain

hist(rain)

plot(qgamma((1/228)*(1:227),shape=.471,scale=1),rain[order(rain)])


#The empirical  cumulative distribution function (chapter 10)
#Put the BEESWAX data set in a convenient place

#beeswax<-read.table("U:/Math 152/BEESWAX.txt",header=T)
beeswax<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/BEESWAX.txt",header=T)


attach(beeswax)
names(beeswax)
#[1] "mp"  "phc"

plot(ecdf(mp))

About 90% of samples have melting point less than 64.2 C and about 12% had melting point below
63.2C. The investigators determined that an addition of 5% microcrystalline wax raises the 
melting point by .85C and that addition of 10% microcrystalline wax
 raises the boiling point by 2.22C.  We see from the plot of the ecdf that the 
10% addition would be easily detectable.



# Chapter 10 Problem 2

l<-runif(16)
g<-function(x) {ecdf(l)(x) -x}
 plot(g,ylim=c(-1,1))
 curve(sqrt(x*(1-x)/16),0,1,add=T)
 curve(-sqrt(x*(1-x)/16),0,1,add=T)
 curve(2*sqrt(x*(1-x)/16),0,1,add=T)
 curve(-2*sqrt(x*(1-x)/16),0,1,add=T)
  curve(-3*sqrt(x*(1-x)/16),0,1,add=T)
 curve(3*sqrt(x*(1-x)/16),0,1,add=T)

l<-runif(16)
g<-function(x) {ecdf(l)(x) -x}
 plot(g,ylim=c(-1,1),add=T)

#chapter 10 problem 8


l<--log(runif(100))

curve(log(1-ecdf(l)(x))+x,0,20,ylim=c(-1.5,1.5))

curve(.1*sqrt(exp(x)-1),0,10,add=T)
curve(.2*sqrt(exp(x)-1),0,10,add=T)
curve(.3*sqrt(exp(x)-1),0,10,add=T)
curve(-.1*sqrt(exp(x)-1),0,10,add=T)
curve(-.2*sqrt(exp(x)-1),0,10,add=T)
curve(-.3*sqrt(exp(x)-1),0,10,add=T)

l<--log(runif(100))
curve(log(1-ecdf(l)(x))+x,0,10,add=T)


#Guinea Pig data

#gpigs<-read.table("U:/Math 152/GUINEPIG.txt",header=T)
gpigs<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/GUINEPIG.txt",header=T,na.strings="*")



attach(gpigs)
names(gpigs)
gpigs

f<-function(list){list[!is.na(list)]}
C<-f(C)
I<-f(I)
II<-f(II)
III<-f(III)
IV<-f(IV)
V<-f(V)

Survival functions:

 curve(1-ecdf(C)(x),0,750)
 curve(1-ecdf(I)(x),0,750,lty=2,add=T)
 curve(1-ecdf(II)(x),0,750,lty=3,add=T)
 curve(1-ecdf(III)(x),0,750,lty=4,add=T)
 curve(1-ecdf(IV)(x),0,750,lty=5,add=T)
 curve(1-ecdf(V)(x),0,750,lty=6,add=T)

Hazard functions: 

curve(log(1-ecdf(C)(x)),0,750)
 curve(log(1-ecdf(I)(x)),0,750,lty=2,add=T)
 curve(log(1-ecdf(II)(x)),0,750,lty=3,add=T)
 curve(log(1-ecdf(III)(x)),0,750,lty=4,add=T)
 curve(log(1-ecdf(IV)(x)),0,750,lty=5,add=T)
 curve(log(1-ecdf(V)(x)),0,750,lty=6,add=T)

A Q-Q plot:

op<-par(no.readonly=TRUE)

par(mfcol=c(2,1))


plot(V[order(V)],III[order(III)])
abline(0,1)

curve(1-ecdf(III)(x),0,750,lty=4)
curve(1-ecdf(V)(x),0,750,lty=6,add=T)

par(op)
For the weakest 60%, the treatment effect is multiplicative. 
Lower dose animals live about twice as long. For the strongest 40%, the treatment effect is 
roughly additive. The lower dose animals tend to live about 100 days longer.

#Heat of sublimation of platinum, pg 361

#platinum<-read.table("U:/Math 152/PLATINUM.txt",header=T)
platinum<-read.table("http://math.cmc.edu/moneill/Math152/Handouts/PLATINUM.txt",header=T)

attach(platinum)
names(platinum)

#A naive confidence interval for the sample mean:

c(mean(hsub) -(qnorm(.975))*sd(hsub)/sqrt(length(hsub)),mean(hsub) +(qnorm(.975))*sd(hsub)/sqrt(length(hsub)))

summary(hsub)

plot(1:26,hsub)

#trimmed mean:

mean(hsub,trim=.2)

#M-estimates:

library(MASS)

?huber

huber(hsub)

#Bootstrapping for the variability of location estimates:

op<-par(no.readonly=TRUE)
par(mfcol=c(4,1),pin=c(1.8,.9))

d<-replicate(1000,mean(sample(hsub,replace=T)))
hist(d,xlim=c(134,140))

d<-replicate(1000,mean(sample(hsub,replace=T),trim=.2))
hist(d,xlim=c(134,140))

d<-replicate(1000,median(sample(hsub,replace=T)))
hist(d,xlim=c(134,140))

d<-replicate(1000,as.numeric(huber(sample(hsub,replace=T))[1]))
hist(d,xlim=c(134,140))


par(op)

#Bootstrapping for confidence intervals for the trimmed mean and median

d<-replicate(1000,mean(sample(hsub,replace=T),trim=.2))
thetahat<-mean(hsub,trim=.2)

deltahigh<-d[order(d)][950]-thetahat
deltalow<-d[order(d)][50]-thetahat

c(thetahat-deltahigh,thetahat-deltalow)

d<-replicate(1000,median(sample(hsub,replace=T)))
thetahat<-median(hsub)

deltahigh<-d[order(d)][950]-thetahat
deltalow<-d[order(d)][50]-thetahat

c(thetahat-deltahigh,thetahat-deltalow)

#Robust measures of dispersion

?IQR

?mad

mad(hsub)

median(abs(hsub-median(hsub)))*1.4826

#----------------------------------------------------------------------
#Boxplots and t-tests

#intrinsic vs. extrinsic motivation in creative writers

#haiku<-read.csv("U:/Math 152/haiku.txt",header=T)
#haiku<-read.csv("C:/Documents and Settings/Mike/Desktop/Math 152/haiku.txt",header=T)
haiku<-read.csv("http://math.cmc.edu/moneill/Math152/Handouts/haiku.txt",header=T)

attach(haiku)
names(haiku)

score.int<-SCORE[TREATMENT=="INTRINSIC"]

score.ext<-SCORE[TREATMENT=="EXTRINSIC"]

boxplot(SCORE~TREATMENT)

?boxplot


#--------------------------------------------



# from example A (and B) section 11.2

A<-c(79.98,80.04,80.02,80.04,80.03,80.03,80.04,79.97,80.05,80.03,80.02,80.00,80.02)
 B<-c(80.02,79.94,79.98,79.97,79.97,80.03,79.95,79.97)

boxplot(A,B)

#exercise 39 pg 414
(remove the "&"'s from the data set then read it in as a matrix)

#winds<-read.table("U:/Math 152/WINDS.txt",header=T)
#winds<-matrix(scan("C:/Documents and Settings/Mike/Desktop/Math 152/WINDS.txt",n=21*37),21,37,byrow=TRUE)
winds<-matrix(scan("http://math.cmc.edu/moneill/Math152/Handouts/WINDS.txt",n=21*37),21,37,byrow=TRUE)
boxplot(data.frame(winds))
boxplot(data.frame(t(winds)))