STAT W4702 - UNI-MM3557 - Assignment 10

From Introduction to Statistical Learning using R

5.4.1

#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#

5.4.2 a

#
#
#
#
#
#
#
#

#
#
#
#
#
#
#
#

#
#
#
#
#
#
#
#

#
#
#
#
#
#
#
#


e
#
#
#
#
#
#
#
#

f
#
#
#
#
#
#
#
#

pr = function(n) return(1 - (1 - 1/n)^n)
x = 1:1e+05
plot(x, pr(x))

5.4.6

library(ISLR)
summary(Default)

##  default    student       balance           income     
##  No :9667   No :7056   Min.   :   0.0   Min.   :  772  
##  Yes: 333   Yes:2944   1st Qu.: 481.7   1st Qu.:21340  
##                        Median : 823.6   Median :34553  
##                        Mean   : 835.4   Mean   :33517  
##                        3rd Qu.:1166.3   3rd Qu.:43808  
##                        Max.   :2654.3   Max.   :73554

set.seed(1)

glm.fit = glm(default ~ income + balance, data = Default, family = binomial)
summary(glm.fit)

## 
## Call:
## glm(formula = default ~ income + balance, family = binomial, 
##     data = Default)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.4725  -0.1444  -0.0574  -0.0211   3.7245  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept) -1.154e+01  4.348e-01 -26.545  < 2e-16 ***
## income       2.081e-05  4.985e-06   4.174 2.99e-05 ***
## balance      5.647e-03  2.274e-04  24.836  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2920.6  on 9999  degrees of freedom
## Residual deviance: 1579.0  on 9997  degrees of freedom
## AIC: 1585
## 
## Number of Fisher Scoring iterations: 8

boot.fn = function(data, index) return(coef(glm(default ~ income + balance, data = data, family = binomial, subset = index)))

library(boot)
boot(Default, boot.fn, 50)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Default, statistic = boot.fn, R = 50)
## 
## 
## Bootstrap Statistics :
##          original        bias     std. error
## t1* -1.154047e+01  1.181200e-01 4.202402e-01
## t2*  2.080898e-05 -5.466926e-08 4.542214e-06
## t3*  5.647103e-03 -6.974834e-05 2.282819e-04

The results are approximatly the same in both cases.

5.4.9

library(MASS)
summary(Boston)

##       crim                zn             indus            chas        
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Min.   :0.00000  
##  1st Qu.: 0.08204   1st Qu.:  0.00   1st Qu.: 5.19   1st Qu.:0.00000  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Median :0.00000  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14   Mean   :0.06917  
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10   3rd Qu.:0.00000  
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74   Max.   :1.00000  
##       nox               rm             age              dis        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##       rad              tax           ptratio          black       
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##      lstat            medv      
##  Min.   : 1.73   Min.   : 5.00  
##  1st Qu.: 6.95   1st Qu.:17.02  
##  Median :11.36   Median :21.20  
##  Mean   :12.65   Mean   :22.53  
##  3rd Qu.:16.95   3rd Qu.:25.00  
##  Max.   :37.97   Max.   :50.00

set.seed(1)
attach(Boston)

medv.mean = mean(medv)
medv.mean

## [1] 22.53281

medv.err = sd(medv)/sqrt(length(medv))
medv.err

## [1] 0.4088611

The results are approximatly the sam eas in the case of b.

boot.fn = function(data, index) return(mean(data[index]))
library(boot)
bstrap = boot(medv, boot.fn, 1000)
bstrap

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original      bias    std. error
## t1* 22.53281 0.008517589   0.4119374

t.test(medv)

## 
##  One Sample t-test
## 
## data:  medv
## t = 55.111, df = 505, p-value < 2.2e-16
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  21.72953 23.33608
## sample estimates:
## mean of x 
##  22.53281

c(bstrap$t0 - 2 * 0.4119, bstrap$t0 + 2 * 0.4119)

## [1] 21.70901 23.35661

medv.med = median(medv)
medv.med

## [1] 21.2

boot.fn = function(data, index) return(median(data[index]))
boot(medv, boot.fn, 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*     21.2 -0.0098   0.3874004

we observe a smaller standard error as comapred to the median.

medv.tenth = quantile(medv, c(0.1))
medv.tenth

##   10% 
## 12.75

boot.fn = function(data, index) return(quantile(data[index], c(0.1)))
boot(medv, boot.fn, 1000)

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = medv, statistic = boot.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##     original  bias    std. error
## t1*    12.75 0.00515   0.5113487

From Introduction to Mathematical Statistics:

4.9.1

The bootstrap sample is i.i.d as the draw is with replacement ```{}

c # # # # # # # #

d # # # # # # # #

4.9.3

b. The updated code: percentcibootMedian<-function(x,b,alpha){ # x is a vector containing the original sample. # b is the desired number of bootstraps. # alpha: (1 - alpha) is the confidence coefficient. # # theta is the # lower is the # upper is the # thetastar is # theta<-median(x) thetastar<-rep(0,b) n<-length(x) for(i in 1:b){xstar<-sample(x,n,replace=T) thetastar[i]<-mean(xstar) } thetastar<-sort(thetastar) pick<-round((alpha/2)*(b+1)) lower<-thetastar[pick] upper<-thetastar[b-pick+1] list(theta=theta,lower=lower,upper=upper,thetastar=thetastar) #list(theta=theta,lower=lower,upper=upper) } ```

4.9.5

The updated code:

newcode<-function(x,b,alpha){
  # x is a vector containing the original sample.
# b is the desired number of bootstraps.
# alpha: (1 - alpha) is the confidence coefficient. #
# theta is the
# lower is the
# upper is the
# thetastar is
#
theta<-mean(x) 
stan<-var(x)^.5
teestar<-rep(0,b)
n<-length(x)
for(i in 1:b){xstar<-sample(x,n,replace=T)
theestar[i]=(mean(xstar)-theta)/(var(xstar)^.5/sqrt(n))
} 
teestar<-sort(teestar)
pick<-round((alpha/2)*(b+1))
lower0<-teestar[pick]
upper0<-teestar[b-pick+1] 
lower<-theta-upper0*(stan/sqrt(n))
upper<-theta-lower0*(stan/sqrt(n)) 
list(theta=theta,lower=lower,upper=upper,teestar=teestar) 
#list(theta=theta,lower=lower,upper=upper)
}

4.9.10

#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#


4.9.13
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#
#

STAT W4702 - UNI-MM3557 - Assignment 10

Mayank Misra

December 1, 2015