bootstrap-and-jackknife Solutions 2022

bootstrap-and-jackknife–1-.R
VIP
2020-04-16# Question No 2a.
# A btstrap sampling estimates for 100 samples with replacement
btsampl <- function(x,sampNo=100){
x = as.matrix(x)
nx = nrow(x)
btsamp = replicate(sampNo,x[sample.int(nx,replace=TRUE),])}
btse <- function(btsamp,fun,...)
{if(is.matrix(btsamp)){
theta = apply(btsamp,2,fun,...)
} else {
theta = apply(btsamp,3,fun,...)}
if(is.matrix(theta)){
return(list(theta=theta,cov=cov(t(theta))))}
else{
return(list(theta=theta,se=sd(theta)))}}
btbias <- function(btse,theta,...){
if(is.matrix(btse$theta)){
return(apply(btse$theta,1,mean) - theta)}
else{return(mean(btse$theta) - theta)}}
# From the data below, the estimate for variance using is given by:
# X~N(20,mean=1,var=1)( Helwig, 2017
)
x = rnorm(20,1,1)
btsamp = btsampl(x)
btse = btse(btsamp,mean)
var(x)
## [1] 1.431905
btse$se
## [1] 0.2260776
mean(btsamp)
## [1] 1.003066
sd(btsamp)
## [1] 1.157066
hist(btse$theta)

#----------------------------------------------------------
# Calculation for the jkknife set up
jksamp <- function(x){
nx = length(x)
jsamp = matrix(0,nx-1,nx)
for(j in 1:nx) jsamp[,j] = x[-j]
jsamp}
jkse <- function(jsamp,fun,...){
nx = ncol(jsamp)
theta = apply(jsamp,2,fun,...)
se = sqrt( ((nx-1)/nx)*sum( (theta-mean(theta))^2 ) )
list(theta=theta,se=se)}
# Calculation for the jkknife estimates for variance
# X~N(20,mean=1,var=1)( Chernick, 2012
)
x = rnorm(20,1,1)
jsamp = jksamp(x)
jse = jkse(jsamp,mean)
var(x)
## [1] 1.664876
sd(x)/sqrt(20)
## [1] 0.28852
jse$se
## [1] 0.28852

mean(jsamp)
## [1] 0.9336294
sd(jsamp)
## [1] 1.259288
hist(jse$theta)
#=====================================================================
=
#Question 2b.
#for btstrap calculations
btsampl <- function(x,sampNo=100){
x = as.matrix(x)
nx = nrow(x)
btsamp = replicate(sampNo,x[sample.int(nx,replace=TRUE),])}
btse <- function(btsamp,fun,...)
{if(is.matrix(btsamp)){
theta = apply(btsamp,2,fun,...)
} else {
theta = apply(btsamp,3,fun,...)}
if(is.matrix(theta)){
return(list(theta=theta,cov=cov(t(theta))))}
else{
return(list(theta=theta,se=sd(theta)))}}
btbias <- function(btse,theta,...){