# Problem_10_9_26.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
par(mfcol=c(2,2))
hist(x,main="Platinum")
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)

x=x.rhodium
par(mfcol=c(2,2))

hist(x,main="Rhodium")
stem(x)
##
##   The decimal point is at the |
##
##   126 | 4
##   127 |
##   128 |
##   129 |
##   130 |
##   131 | 111112234569
##   132 | 123456677899
##   133 | 000333455558
##   134 | 12
##   135 | 7
boxplot(x)
plot(x)

x=x.iridium
par(mfcol=c(2,2))

hist(x,main="Iridium")
stem(x)
##
##   The decimal point is 1 digit(s) to the right of the |
##
##   13 | 7
##   14 |
##   14 | 5
##   15 | 2
##   15 | 999
##   16 | 00000000000000001113
##   16 |
##   17 | 4
boxplot(x)
plot(x)
plot(c(10:length(x)), x[10:length(x)], main="Iridium")

# (e) For Platinum, observations 8,9,10, and 14,5 seem much higher than the rest.
#     They do not seem iid

# (e) For Rhodium,  the first 2 observations are very variable (low then high).
#  and the observations after about number 25 seem different from those before.
#  These data do not seem iid either.

# (e) For Iridium,  the first 4 observations steadily increase and observaton 8 seems much higher
# than the rest.  These data do not appear iid. Also, plotting the observations from 10 on,
# there seems to be time series dependence -- closer observations in the sequence have more similar
# values.

# (f) Measures of location
#   For Rhodium
x=x.rhodium
mean(x); mean(x,trim=.1);mean(x,trim=.2);median(x); huber(x,k=1.5)[[1]]
## [1] 132.42
## [1] 132.4781
## [1] 132.5292
## [1] 132.65
## [1] 132.4921
# All the estimates are very close
# The mean is the lowest, it is affected by the lowest value
# The trimmed means exclude the extreme values and are similar to the median.

#   For Iridium
x=x.iridium
mean(x); mean(x,trim=.1);mean(x,trim=.2);median(x); huber(x,k=1.5)[[1]]
## [1] 158.8148
## [1] 159.5478
## [1] 159.8412
## [1] 159.8
## [1] 159.8545
# The mean is much lower than the other estimates
# The mean is the lowest, it is affected by the lowest value
# The trimmed means exclude the extreme values and are similar to the median.

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

print(x.mean.sterr)
## [1] 0.2273369
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] 132.42
## [1] 132.0461 132.7939
stats1.rhodium=c(x.mean, x.mean.ci)

# Second for iridium
x=x.iridium
x.stdev=sqrt(var(x))
x.mean.sterr=x.stdev/sqrt(length(x))

print(x.mean.sterr)
## [1] 1.197917
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] 158.8148
## [1] 156.8444 160.7852
stats1.iridium=c(x.mean, x.mean.ci)

# 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]])
}

#  First for rhodium
set.seed(1)
x=x.rhodium

# Bootstrap analysis of sample 20% trimmed mean:
x.boot.trimmedmean.2= boot(x, fcn.trimmedmean,trim=.2, R=1000)
par(mfcol=c(1,2))
plot(x.boot.trimmedmean.2)
title(paste("Trimmed Mean:  StDev=",as.character(round(sqrt(var(x.boot.trimmedmean.2$t)),digits=4)),sep=""), adj=1) print(x.boot.trimmedmean.2$t0)
## [1] 132.42
print(sqrt(var(x.boot.trimmedmean.2$t))) ## [,1] ## [1,] 0.2214742 stats2.trim20.rhodium<-c(x.boot.trimmedmean.2$t0,(sqrt(var(x.boot.trimmedmean.2$t)))) # boot.ci(x.boot.trimmedmean.2,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.trimmedmean.2, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (132.1, 132.8 ) ## Calculations and Intervals on Original Scale # Bootstrap analysis of sample 10% trimmed mean: x.boot.trimmedmean.1= boot(x, fcn.trimmedmean,trim=.1, R=1000) par(mfcol=c(1,2)) plot(x.boot.trimmedmean.1) title(paste("Trimmed Mean: StDev=",as.character(round(sqrt(var(x.boot.trimmedmean.1$t)),digits=4)),sep=""),
adj=1)

print(x.boot.trimmedmean.1$t0) ## [1] 132.42 print(sqrt(var(x.boot.trimmedmean.1$t)))
##           [,1]
## [1,] 0.2195023
stats2.trim10.rhodium<-c(x.boot.trimmedmean.1$t0,(sqrt(var(x.boot.trimmedmean.1$t))))
#
#
boot.ci(x.boot.trimmedmean.1,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.trimmedmean.1, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (132.1, 132.8 )
## Calculations and Intervals on Original Scale
# 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] 132.65
print(sqrt(var(x.boot.median$t))) ## [,1] ## [1,] 0.2090489 boot.ci(x.boot.median,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.median, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (132.3, 133.0 ) ## Calculations and Intervals on Original Scale stats2.median.rhodium<-c(x.boot.median$t0,(sqrt(var(x.boot.median$t)))) stats1.rhodium ## [1] 132.4200 132.0461 132.7939 # The bootstrap estimate and standard errrors are given by: stats2.trim20.rhodium ## [1] 132.4200000 0.2214742 stats2.trim10.rhodium ## [1] 132.4200000 0.2195023 stats2.median.rhodium ## [1] 132.6500000 0.2090489 # For Rhodium these are all about the same # (k). The approximate 90% confidence intervals are: boot.ci(x.boot.trimmedmean.2,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.trimmedmean.2, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (132.1, 132.8 ) ## Calculations and Intervals on Original Scale boot.ci(x.boot.trimmedmean.1,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.trimmedmean.1, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (132.1, 132.8 ) ## Calculations and Intervals on Original Scale boot.ci(x.boot.median,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.median, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (132.3, 133.0 ) ## Calculations and Intervals on Original Scale # These are all essential equal for rhodium # These intervals are all about the same as that based on part (g) stats1.rhodium ## [1] 132.4200 132.0461 132.7939 # ############## # Second for iridium set.seed(1) x=x.iridium # Bootstrap analysis of sample 20% trimmed mean: x.boot.trimmedmean.2= boot(x, fcn.trimmedmean,trim=.2, R=1000) par(mfcol=c(1,2)) plot(x.boot.trimmedmean.2) title(paste("Trimmed Mean: StDev=",as.character(round(sqrt(var(x.boot.trimmedmean.2$t)),digits=4)),sep=""),
adj=1)

print(x.boot.trimmedmean.2$t0) ## [1] 158.8148 print(sqrt(var(x.boot.trimmedmean.2$t)))
##          [,1]
## [1,] 1.181086
stats2.trim20.iridium<-c(x.boot.trimmedmean.2$t0,(sqrt(var(x.boot.trimmedmean.2$t))))
#
boot.ci(x.boot.trimmedmean.2,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.trimmedmean.2, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (157.0, 160.8 )
## Calculations and Intervals on Original Scale
# Bootstrap analysis of sample 10% trimmed mean:
x.boot.trimmedmean.1= boot(x, fcn.trimmedmean,trim=.1, R=1000)
par(mfcol=c(1,2))
plot(x.boot.trimmedmean.1)
title(paste("Trimmed Mean:  StDev=",as.character(round(sqrt(var(x.boot.trimmedmean.1$t)),digits=4)),sep=""), adj=1) print(x.boot.trimmedmean.1$t0)
## [1] 158.8148
print(sqrt(var(x.boot.trimmedmean.1$t))) ## [,1] ## [1,] 1.185398 stats2.trim10.iridium<-c(x.boot.trimmedmean.1$t0,(sqrt(var(x.boot.trimmedmean.1$t)))) # # boot.ci(x.boot.trimmedmean.1,conf=.90, type="basic") ## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS ## Based on 1000 bootstrap replicates ## ## CALL : ## boot.ci(boot.out = x.boot.trimmedmean.1, conf = 0.9, type = "basic") ## ## Intervals : ## Level Basic ## 90% (157.0, 160.9 ) ## Calculations and Intervals on Original Scale # 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] 159.8 print(sqrt(var(x.boot.median$t)))
##           [,1]
## [1,] 0.2127957
boot.ci(x.boot.median,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.median, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (159.5, 160.1 )
## Calculations and Intervals on Original Scale
stats2.median.iridium<-c(x.boot.median$t0,(sqrt(var(x.boot.median$t))))

stats1.iridium
## [1] 158.8148 156.8444 160.7852
stats2.trim20.iridium
## [1] 158.814815   1.181086
stats2.trim10.iridium
## [1] 158.814815   1.185398
stats2.median.iridium
## [1] 159.8000000   0.2127957
##############
# The bootstrap estimate and standard errrors are given by:
stats2.trim20.iridium
## [1] 158.814815   1.181086
stats2.trim10.iridium
## [1] 158.814815   1.185398
stats2.median.iridium
## [1] 159.8000000   0.2127957
# For iridium the median has much lower st error than the trimmed means

#  (k). The approximate 90% confidence intervals are:

boot.ci(x.boot.trimmedmean.2,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.trimmedmean.2, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (157.0, 160.8 )
## Calculations and Intervals on Original Scale
boot.ci(x.boot.trimmedmean.1,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.trimmedmean.1, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (157.0, 160.9 )
## Calculations and Intervals on Original Scale
boot.ci(x.boot.median,conf=.90, type="basic")
## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
## Based on 1000 bootstrap replicates
##
## CALL :
## boot.ci(boot.out = x.boot.median, conf = 0.9, type = "basic")
##
## Intervals :
## Level      Basic
## 90%   (159.5, 160.1 )
## Calculations and Intervals on Original Scale
# The interval using the median is much smaller than that for the trimmed means
# These intervals are all smaller than that based on part (g)
stats1.iridium
## [1] 158.8148 156.8444 160.7852