Added Bootstrap Material

material/4_thu/bootstrap/CodeSnippets.jl

@@ -0,0 +1,104 @@
+using Distributions
+using Statistics
+using StatsPlots
+
+d = Normal(0.0, 1.0)
+n = 100
+x = rand(d, n)
+θ = mean(x)
+
+#---
+
+σ = std(x)
+est_d = Normal(θ, σ/sqrt(n))
+plot(est_d, legend=false)
+
+ci_bounds = quantile(est_d, [0.025,0.975])
+vline!(ci_bounds)
+
+#---
+
+B = 1000
+est_vector = zeros(B)
+for i in 1:B
+  x_new = rand(d, n)
+  est_vector[i] = mean(x_new)
+end    
+histogram(est_vector, legend=false)
+
+ci_bounds_new = quantile(est_vector, [0.025, 0.975])
+vline!(ci_bounds_new)
+
+#---
+
+est_vector_bs = zeros(B)
+for i in 1:B
+  x_bs = rand(x, n)
+  est_vector_bs[i] = mean(x_bs)
+end
+histogram(est_vector_bs, legend=false)
+
+ci_bounds_bs = quantile(est_vector_bs, [0.025, 0.975])
+vline!(ci_bounds_bs) 
+
+#---
+
+# the following function generates a sample of (strongly) 
+# correlated normal data with mean μ and variance 1
+myrandn = function (mu=0.0, n=1)
+  rho = 0.9
+  res = zeros(n)
+  res[1] = randn(1)[1]
+  if n > 1
+    for i in 2:n
+      res[i] = rho*res[i-1] .+ sqrt(1-rho^2)*randn(1)[1]
+    end
+  end
+  res .= mu .+ res
+  return(res)
+end 
+
+## Assume we obtained the following data
+n = 20
+true_mu = 1.0
+y = myrandn(true_mu, n)
+est_mu = mean(y)
+
+## theoretical standard error
+B = 1000
+est_vector = zeros(B)
+for i in 1:B
+  y_new = myrandn(true_mu, n)
+  est_vector[i] = mean(y_new)
+end    
+std(est_vector)
+
+## classical i.i.d. bootstrap
+est_vector_bs = zeros(B)
+for i in 1:B
+  y_bs = rand(y, n)
+  est_vector_bs[i] = mean(y_bs)
+end
+std(est_vector_bs)
+
+## parametric bootstrap
+est_vector_pbs = zeros(B)
+for i in 1:B
+  y_pbs = myrandn(est_mu, n)
+  est_vector_pbs[i] = mean(y_pbs)
+end    
+std(est_vector_pbs)
+
+# ---
+myrand = function(mu, n) 
+  rho = 0.9
+  res = zeros(n)
+  res[1] = rand(1)[1]
+  if n > 1
+    for i in 2:n
+      res[i] = rho*res[i-1] .+ (1-rho)*rand(1)[1]
+    end
+  end
+  res .= mu - 0.5 .+ res
+  return(res)
+end