2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 1/13
5291 HW1
Mengjia Huang, UNI:mh3781
2018/9/13
Data that remove the column of “dose”
library(datasets)
data(ToothGrowth)
head(ToothGrowth)
## len supp dose
## 1 4.2 VC 0.5
## 2 11.5 VC 0.5
## 3 7.3 VC 0.5
## 4 5.8 VC 0.5
## 5 6.4 VC 0.5
## 6 10.0 VC 0.5
data=subset(ToothGrowth,select = c("len","supp"))
i)Perform. separate EDA for each supplement , (i.e., OJ /VC).,and compute the Sample kurtosis, Sample median,Sample Inter
Quartile Range (IQR), and Sample IQR for len in each group.
#make sure whether there is missing data
sum(is.na(data))
## [1] 0
VC=data[data$supp=="VC",]
OJ=data[data$supp=="OJ",]
EDA for group OJ
summary(OJ$len)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 8.20 15.52 22.70 20.66 25.72 30.90
stem(OJ$len)
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 2/13
##
## The decimal point is at the |
##
## 8 | 2477
## 10 | 0
## 12 |
## 14 | 552
## 16 | 56
## 18 | 7
## 20 | 025
## 22 | 4036
## 24 | 58258
## 26 | 44433
## 28 | 4
## 30 | 9
var(OJ$len)
## [1] 43.63344
sd(OJ$len)
## [1] 6.605561
IQR(OJ$len)
## [1] 10.2
boxplot(OJ$len, main="boxplot of len in group OJ ")
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 3/13
hist(OJ$len,main="histogram of len in group OJ")
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 4/13
qqnorm(OJ$len,main="qqplot of len in group OJ")
qqline(OJ$len)
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 5/13
Sample kurtosis for group OJ
library(e1071)
## Warning: package 'e1071' was built under R version 3.3.2
library(bootstrap)
## Warning: package 'bootstrap' was built under R version 3.3.2
kurtosis(OJ$len)
## [1] -1.030564
Sample median for group OJ
median(OJ$len)
## [1] 22.7
Sample Inter Quartile Range (IQR) for group OJ
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 6/13
IQR(OJ$len)
## [1] 10.2
EDA for group VC
summary(VC$len)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.20 11.20 16.50 16.96 23.10 33.90
stem(VC$len)
##
## The decimal point is 1 digit(s) to the right of the |
##
## 0 | 4
## 0 | 56677
## 1 | 01124
## 1 | 556777799
## 2 | 2334
## 2 | 667
## 3 | 034
var(VC$len)
## [1] 68.32723
sd(VC$len)
## [1] 8.266029
IQR(VC$len)
## [1] 11.9
boxplot(VC$len, main="boxplot of len in group VC ")
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 7/13
hist(VC$len,main="histogram of len in group VC")
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 8/13
qqnorm(VC$len,main="qqplot of len in group VC")
qqline(VC$len)
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 9/13
sample kurtosis for group VC
kurtosis(VC$len)
## [1] -0.9274193
Sample median for group VC
median(VC$len)
## [1] 16.5
Sample Inter Quartile Range (IQR) for group VC
IQR(VC$len)
## [1] 11.9
ii)For each of the estimates computed in (i) above, determinethe bias and variance using methods of Jackknife and
Bootstrap.
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 10/13
Jackknife for group OJ
library(bootstrap)
jackknife(OJ$len,kurtosis)$jack.bias #bias of kurtosis
## [1] -0.01569017
(jackknife(OJ$len,kurtosis)$jack.se)^2 #variance of kurtosis
## [1] 0.2427639
jackknife(OJ$len,median)$jack.bias #bias of median
## [1] 0
(jackknife(OJ$len,median)$jack.se)^2 #variacne of kurtosis
## [1] 2.61
jackknife(OJ$len,IQR)$jack.bias #bias of IQR
## [1] -0.7733333
(jackknife(OJ$len,IQR)$jack.se)^2 #variance of IQR
## [1] 8.486044
Jackknife for group VC
jackknife(VC$len,kurtosis)$jack.bias #bias of kurtosis
## [1] -0.1174358
(jackknife(VC$len,kurtosis)$jack.se)^2 #variance of kurtosis
## [1] 0.1227792
jackknife(VC$len,median)$jack.bias #bias of median
## [1] 0
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 11/13
(jackknife(VC$len,median)$jack.se)^2 #variacne of kurtosis
## [1] 0
jackknife(VC$len,IQR)$jack.bias #bias of IQR
## [1] -0.3866667
(jackknife(VC$len,IQR)$jack.se)^2 #variance of IQR
## [1] 3.629511
Bootstrap for group OJ
set.seed(0)
B=1000
n=length(OJ$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(OJ$len,n,replace = TRUE)
est[i]=kurtosis(resample)
}
B_Kurt_bias=mean(est)-kurtosis(OJ$len) #bias of Kurtosis
B_Kurt_bias
## [1] 0.1212991
B_Kurt_var=var(est)#var of Kurtosis
B_Kurt_var
## [1] 0.2674562
set.seed(0)
B=1000
n=length(OJ$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(OJ$len,n,replace = TRUE)
est[i]=median(resample)
}
B_Med_bias=mean(est)-median(OJ$len) #bias of Median
B_Med_bias
## [1] -0.36905
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 12/13
B_Med_var=var(est)#var of Median
B_Med_var
## [1] 3.033728
set.seed(0)
B=1000
n=length(OJ$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(OJ$len,n,replace = TRUE)
est[i]=IQR(resample)
}
B_IQR_bias=mean(est)-IQR(OJ$len) #bias of IQR
B_IQR_bias
## [1] -0.5843
B_IQR_var=var(est)#var of IQR
B_IQR_var
## [1] 6.602072
Bootstrap for group VC
set.seed(0)
B=1000
n=length(VC$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(VC$len,n,replace = TRUE)
est[i]=kurtosis(resample)
}
B_Kurt_bias1=mean(est)-kurtosis(VC$len) #bias of Kurtosis
B_Kurt_bias1
## [1] 0.04572266
B_Kurt_var1=var(est)#var of Kurtosis
B_Kurt_var1
## [1] 0.1358673
2018/9/13 5291 HW1
file://localhost/Users/minta/Desktop/5291/5291_HW1.html 13/13
set.seed(0)
B=1000
n=length(VC$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(VC$len,n,replace = TRUE)
est[i]=median(resample)
}
B_Med_bias1=mean(est)-median(VC$len) #bias of Median
B_Med_bias1
## [1] -0.13145
B_Med_var1=var(est)#var of Median
B_Med_var1
## [1] 3.010559
set.seed(0)
B=1000
n=length(VC$len)
est<-vector(length=B)
for (i in 1:B){
resample=sample(VC$len,n,replace = TRUE)
est[i]=IQR(resample)
}
B_IQR_bias1=mean(est)-IQR(VC$len) #bias of IQR
B_IQR_bias1
## [1] -0.186775
B_IQR_var1=var(est)#var of IQR
B_IQR_var1
## [1] 7.747735