From d6b3f6c2deab482c808ec77f1b4988ab6f6ef84d Mon Sep 17 00:00:00 2001
From: Joe Fowler <drjoefowler@gmail.com>
Date: Tue, 7 Feb 2017 14:33:52 -0700
Subject: [PATCH] Unify README on bootstrapci and bootstrapse()

Bootstrap methods to extract a 1-alpha confidence interval and the
standard error of an estimator.
---
 README.md          | 16 ++++++++-----
 src/RobustStats.jl |  2 +-
 src/functions.jl   | 59 +++++++++++++++++++---------------------------
 test/runtests.jl   | 14 ++++++++---
 4 files changed, 46 insertions(+), 45 deletions(-)

diff --git a/README.md b/README.md
index 2433cc4..2f06148 100644
--- a/README.md
+++ b/README.md
@@ -211,13 +211,17 @@ Compute one-step M-estimator of location using Huber's ψ. The default bending c
 
 
 
-#### `bootstrapci`
-Compute a bootstrap, (1-α) confidence interval for the measure of location corresponding to the argument `est`. By default, a one-step M-estimator of location based on Huber's ψ is used. The default number of bootstrap samples is `nboot=2000`. `nullvalue` is the target value used when computing a p-value. Default α=0.05.
+#### `bootstrapci`, `bootstrapse`
+Compute a bootstrap, (1-α) confidence interval (`bootstrapci`) or a standard error (`bootstrapse`) for the measure of location corresponding to the argument `est`. By default, the median is used. Default α=0.05.
 
-    julia> bootstrapci(x)
-    Estimate:             1.305811
-    Confidence interval:  0.687723       2.259071
-    p value:             < 10e-16
+    julia> ci = bootstrapci(x, est=onestep, nullvalue=0.6)
+     Estimate:             1.305811
+     Confidence interval:  0.687723       2.259071
+     p value:              0.026000
+
+
+    julia> se = bootstrapse(x, est=onestep)
+    0.41956761772722817
 
 
 
diff --git a/src/RobustStats.jl b/src/RobustStats.jl
index ecf529a..79dbae7 100644
--- a/src/RobustStats.jl
+++ b/src/RobustStats.jl
@@ -32,6 +32,7 @@ export
     hpsi,
     onestep,
     bootstrapci,
+    bootstrapse,
     mom,
     momci,
     contam_randn,
@@ -40,7 +41,6 @@ export
     yuend,
 
     trimcibt,
-    bootse,
 
     bisquareWM,
     huberWM,
diff --git a/src/functions.jl b/src/functions.jl
index 2ed1627..e74262e 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -550,7 +550,7 @@ function onestep{S <: Real}(x::AbstractArray{S}, bend::Real=1.28)
     return MED + MAD*A/B
 end
 
-"""`bootstrapci(x; est=onestep, alpha=0.05, nboot=2000)`
+"""`bootstrapci(x; est=onestep, alpha=0.05, nboot=2000, nullvalue=NaN)`
 
 Compute a (1-α) confidence interval for the location-estimator function `est`
 using a bootstrap calculation. The default estimator is `onestep`. If `nullvalue` is
@@ -565,13 +565,9 @@ function bootstrapci{S <: Real}(x::AbstractArray{S}; est::Function=onestep,
     end
     const n = length(x)
     bvec = zeros(nboot)
-    temp = zeros(n)
-    randid=rand(1:n, n, nboot)
     for i = 1:nboot
-        for j = 1:n
-            temp[j] = x[randid[j,i]]
-        end
-        bvec[i]=est(temp)
+        randid=rand(1:n, n)
+        bvec[i]=est(x[randid])
     end
     low::Int = round((alpha/2)*nboot) + 1
     up = nboot-low + 1
@@ -590,6 +586,27 @@ function bootstrapci{S <: Real}(x::AbstractArray{S}; est::Function=onestep,
     output
 end
 
+"""`bootstrapse(x; est=median, alpha=0.05, nboot=2000)`
+
+Compute the standard error of the location-estimator function `est`
+using a bootstrap calculation. The default estimator is `median`.
+"""
+function bootstrapse{S <: Real}(x::AbstractArray{S};
+        nboot::Integer=1000, est::Function=median, seed=2)
+    if isa(seed, Int)
+        srand(seed)
+    elseif seed
+        srand(2)
+    end
+    const n = length(x)
+    bvec = zeros(nboot)
+    for i = 1:nboot
+        randid=rand(1:n, n)
+        bvec[i]=est(x[randid])
+    end
+    std(bvec)
+end
+
 
 """`mom(x; bend=2.24)`
 
@@ -731,34 +748,6 @@ function trimcibt{S <: Real}(x::AbstractArray{S}; tr::Real=0.2, alpha::Real=0.05
     end
 end
 
-#   Compute bootstrap estimate of the standard error of the
-#   estimator est
-#   The default number of bootstrap samples is nboot=1000
-
-function bootse{S <: Real}(x::AbstractArray{S}; nboot::Integer=1000, est::Function=median, seed=2)
-   if isa(seed, Bool)
-        if seed
-            srand(2)
-        end
-    else
-        srand(seed)
-    end
-    n=length(x)
-    temp=zeros(n)
-    bvec=zeros(nboot)
-    randid=rand(1:n, n*nboot)
-    for i=1:(nboot*n)
-        if (i%n)!=0
-            temp[i%n]=x[randid[i]]
-        else
-            temp[n]=x[randid[i]]
-            bvec[div(i, n)]=est(temp)
-
-        end
-    end
-    return std(bvec)
-end
-
 
 """`procb(x, y; seed=2)`
 
diff --git a/test/runtests.jl b/test/runtests.jl
index ec77432..caafe80 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -60,9 +60,17 @@ mci = momci(x, seed=2, nboot=2000, nullvalue=0.6)
 @test mci.ci[2] ≈ 2.120978892222222
 @test mci.p ≈ 0.11099999999999999
 
-ci = sint(x).ci
-@test ci[1] ≈ 0.5474825849287168
-@test ci[2] ≈ 2.3752324887983707
+ci = bootstrapci(x, est=onestep, nullvalue=0.6)
+@test ci.estimate ≈ 1.3058109021286803
+@test ci.ci[1] ≈ 0.6877225172104009
+@test ci.ci[2] ≈ 2.259070870017077
+@test ci.p ≈ 0.026000000000000023
+se = bootstrapse(x, est=onestep)
+@test se ≈ 0.41956761772722817
+
+s = sint(x)
+@test s.ci[1] ≈ 0.5474825849287168
+@test s.ci[2] ≈ 2.3752324887983707
 s = sint(x, 0.6)
 @test s.ci[1] ≈ 0.5474825849287168
 @test s.ci[2] ≈ 2.3752324887983707