Skip to content

Instantly share code, notes, and snippets.

@sieste
Last active March 10, 2020 14:29
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save sieste/9740144b77a671f24dfa815fba0100cf to your computer and use it in GitHub Desktop.
Save sieste/9740144b77a671f24dfa815fba0100cf to your computer and use it in GitHub Desktop.
# Example 3.12
x = c(0.02, 0.19, 0.27, 0.28, 0.32, 0.50, 1.08, 1.42, 2.23)
# Example 3.13
x = c(-0.63, 0.18, -0.84, 1.60, 0.33, -0.82, 0.49, 0.74, 0.58)
y = c(-0.31, 1.51, 0.39, -0.62, -2.21, 1.12, -0.04, -0.02, 0.94)
# Example 3.16
x = c(0.38, 0.59, 0.06, 0.08, 0.22, 1.45, 0.61, 0.27, 0.48, 0.07)
B = 1000
alpha = 0.05 # set confidence level
n = length(x)
theta_hat = 1 / mean(x)
t_star = numeric(B)
for (b in 1:B) {
x_star = sample(x, replace=TRUE) # resample X
t_star[b] = 1 / mean(x_star) - theta_hat # calculate T*
}
bas = theta_hat - quantile(t_star, c(1-alpha, alpha)) # basic interval
mle = theta_hat - qnorm(c(1-alpha, alpha)) * theta_hat/sqrt(n) # normal
rbind(mle=mle, basic=bas)
# Example 3.17 (cont)
s = theta_hat / sqrt(n) # estimate standard error
for (b in 1:B) {
x_star = sample(x, replace=TRUE) # resample
theta_star = 1 / mean(x_star) # estimate theta
s_star = theta_star / sqrt(n) # estimate standard error
t_star[b] = (theta_star - theta_hat) / s_star # calculate T*
}
stu = theta_hat - quantile(t_star, c(1-alpha, alpha)) * s # interval
rbind(mle=mle, basic=bas, studentised=stu)
# Example 3.18 (cont)
theta_star = numeric(B) # vector to store estimates
for (b in 1:B) {
x_star = sample(x, replace=TRUE) # resample
theta_star[b] = 1 / mean(x_star) # estimate theta
}
per = quantile(theta_star, c(alpha, 1-alpha)) # percentile interval
rbind(mle=mle, basic=bas, studentised=stu, percentile=per)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment