-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMath312_project2.Rmd
110 lines (99 loc) · 2.46 KB
/
Math312_project2.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
---
title: "312_Project_2"
author: "Kenneth Cho, Thomas Liu, and Andrew MacCausland"
date: "11/21/2021"
output: pdf_document
---
```{r}
library('boot')
```
```{r}
# this function will return the ci object for 1 sample
get_ci <- function(){
# generate 20 samples from standard normal
stdnormal <- rnorm(20)
# calculate the ci object
boot.obj <- boot(stdnormal, R=2000, statistic=function(x, i){mean(x[i])})
}
```
```{r}
set.seed(1331)
norm_miss <- 0
basic_miss <- 0
perc_miss <- 0
# calculate coverage rate, run 1000 times
for (i in 1:1000){
boot.obj = get_ci()
intervals <- boot.ci(boot.obj, type = c("norm", "basic", "perc"))
norm_min <- intervals$normal[2]
norm_max <- intervals$normal[3]
basic_min <- intervals$basic[4]
basic_max <- intervals$basic[5]
perc_min <- intervals$percent[4]
perc_max <- intervals$percent[5]
# get misses
if(0 < norm_min | 0 > norm_max){
norm_miss <- norm_miss + 1
}
if(0 < basic_min | 0 > basic_max){
basic_miss <- basic_miss + 1
}
if(0 < perc_min | 0 > perc_max){
perc_miss <- perc_miss + 1
}
}
cat("coverage rate for standard normal ci: ", 100 - norm_miss / 10, "%\n")
cat("coverage rate for basic ci: ", 100 - basic_miss / 10, "%\n")
cat("coverage rate for percentile ci: ", 100 - perc_miss / 10, "%\n")
```
```{r}
# modify the loop above a bit to get the proportions of times that the
# confidence intervals miss on the left, and the proportion of times
# that the confidence intervals miss on the right
norm_miss_L <- 0
basic_miss_L <- 0
perc_miss_L <- 0
norm_miss_R <- 0
basic_miss_R <- 0
perc_miss_R <- 0
set.seed(1331)
# calculate coverage rate, run 1000 times
for (i in 1:1000){
boot.obj = get_ci()
intervals <- boot.ci(boot.obj, type = c("norm", "basic", "perc"))
norm_min <- intervals$normal[2]
norm_max <- intervals$normal[3]
basic_min <- intervals$basic[4]
basic_max <- intervals$basic[5]
perc_min <- intervals$percent[4]
perc_max <- intervals$percent[5]
# get misses
if(0 < norm_min){
norm_miss_L <- norm_miss_L + 1
}
if(0 < basic_min){
basic_miss_L <- basic_miss_L + 1
}
if(0 < perc_min){
perc_miss_L <- perc_miss_L + 1
}
if(0 > norm_max){
norm_miss_R <- norm_miss_R + 1
}
if(0 > basic_max){
basic_miss_R <- basic_miss_R + 1
}
if(0 > perc_max){
perc_miss_R <- perc_miss_R + 1
}
}
print("normal")
print(norm_miss_L)
print(norm_miss_R)
print("basic")
print(basic_miss_L)
print(basic_miss_R)
print("perc")
print(perc_miss_L)
print(perc_miss_R)
```