# RProject8_3_bootstrap_location.r
#   Problem 10.26 of Rice
#   Hampson and Walker data on heats of sublimation of
#       platinum,iridium, and rhodium

# To install these packages, uncomment the next two lines
#install.packages("MASS")
#install.packages("boot")

library(MASS)

## Warning: package 'MASS' was built under R version 3.1.3

library(boot)

## Warning: package 'boot' was built under R version 3.1.3

x.platinum=scan(file="Rice 3e Datasets/ASCII Comma/Chapter 10/platinum.txt")
x.iridium=scan(file="Rice 3e Datasets/ASCII Comma/Chapter 10/iridium.txt")
x.rhodium=scan(file="Rice 3e Datasets/ASCII Comma/Chapter 10/rhodium.txt")

# Parts (a)-(d)
x=x.platinum
hist(x)

stem(x)

## 
##   The decimal point is at the |
## 
##   132 | 7
##   134 | 134578889900224488
##   136 | 36
##   138 | 
##   140 | 2
##   142 | 3
##   144 | 
##   146 | 58
##   148 | 8

boxplot(x)

plot(x)

# (e)
# Measurements 8,9,10 are all very high as are measurements 14 and 15
# The data do not appear indepenent

# (f) Measures of location
mean(x); mean(x,trim=.1);mean(x,trim=.2);median(x); huber(x,k=1.5)[[1]]

## [1] 137.0462

## [1] 136.3091

## [1] 135.2875

## [1] 135.1

## [1] 135.3841

# (g) Standard error of the sample mean and approximate 90 percent conf. interval
x.stdev=sqrt(var(x))
x.mean.sterr=x.stdev/sqrt(length(x))

print(x.mean.sterr)

## [1] 0.8724542

alpha=0.10
z.upperalphahalf=qnorm(1-alpha/2)
x.mean=mean(x)
x.mean.ci<-c(-1,1)*z.upperalphahalf* x.mean.sterr + x.mean
print(x.mean); print(x.mean.ci)

## [1] 137.0462

## [1] 135.6111 138.4812

# parts (i), (j), (k).
# Bootstrap confidence intervals of different location measures


fcn.median<- function(x, d) {
  return(median(x[d]))
}
fcn.mean<- function(x, d) {
  return(mean(x[d]))
}

fcn.trimmedmean <- function(x, d, trim=0) {
  return(mean(x[d], trim/length(x)))
}

fcn.huber<-function(x,d){
  x.huber=huber(x[d],k=1.5)
  return(x.huber[[1]])
}

# Bootstrap analysis of sample mean:
x.boot.mean= boot(x, fcn.mean, R=1000)      
plot(x.boot.mean)
print(x.boot.mean$t0)

## [1] 137.0462

print(sqrt(var(x.boot.mean$t)))

##           [,1]
## [1,] 0.8406306

boot.ci(x.boot.mean,conf=.95, type="basic")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = x.boot.mean, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (135.2, 138.5 )  
## Calculations and Intervals on Original Scale

title(paste("Sample Mean:  StDev=",as.character(round(sqrt(var(x.boot.mean$t)),digits=4)),sep=""),
      adj=1)

# Bootstrap analysis of sample median:
x.boot.median= boot(x, fcn.median, R=1000)      
plot(x.boot.median)
title(paste("Median:  StDev=",as.character(round(sqrt(var(x.boot.median$t)),digits=4)),sep=""),
      adj=1)

print(x.boot.median$t0)

## [1] 135.1

print(sqrt(var(x.boot.median$t)))

##          [,1]
## [1,] 0.238632

boot.ci(x.boot.median,conf=.95, type="basic")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = x.boot.median, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (134.4, 135.4 )  
## Calculations and Intervals on Original Scale

# Bootstrap analysis of sample trimmed mean:
x.boot.trimmedmean= boot(x, fcn.trimmedmean,trim=.2, R=1000)      
plot(x.boot.trimmedmean)
title(paste("Trimmed Mean:  StDev=",as.character(round(sqrt(var(x.boot.trimmedmean$t)),digits=4)),sep=""),
      adj=1)

print(x.boot.trimmedmean$t0)

## [1] 137.0462

print(sqrt(var(x.boot.trimmedmean$t)))

##          [,1]
## [1,] 0.833879

boot.ci(x.boot.trimmedmean,conf=.95, type="basic")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = x.boot.trimmedmean, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (135.3, 138.6 )  
## Calculations and Intervals on Original Scale

# Bootstrap analysis of Huber M estimate:
x.boot.huber= boot(x, fcn.huber, R=1000)      
plot(x.boot.huber)
title(paste("Huber M Est:  StDev=",as.character(round(sqrt(var(x.boot.huber$t)),digits=4)),sep=""),
      adj=1)

print(x.boot.huber$t0)

## [1] 135.3841

print(sqrt(var(x.boot.huber$t)))

##           [,1]
## [1,] 0.3827326

boot.ci(x.boot.huber,conf=.95, type="basic")

## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
## 
## CALL : 
## boot.ci(boot.out = x.boot.huber, conf = 0.95, type = "basic")
## 
## Intervals : 
## Level      Basic         
## 95%   (134.4, 135.9 )  
## Calculations and Intervals on Original Scale