Merge branch 'master' into two-step-adjust-method

Author: juliangehring

src/pi0-estimators.jl

@@ -68,7 +68,7 @@ end
 function bootstrap_pi0{T<:AbstractFloat,S<:AbstractFloat}(pValues::Vector{T}, lambda::Vector{S} = [0.05:0.05:0.95;], q::S = 0.1)
     #validPValues(pValues)
     n = length(pValues)
-    w = Int[sum(pValues .>= l) for l in lambda] ## TODO: check >= or >
+    w = [sum(pValues .>= l) for l in lambda]  # TODO: check if >= or >
     pi0 = w ./ n ./ (1. - lambda)
     min_pi0 = quantile(pi0, q)
     mse = (w ./ (n.^2 .* (1. - lambda).^2 )) .* (1. - w/n) + (pi0 - min_pi0).^2

src/utils.jl

@@ -10,19 +10,20 @@ end
 
 # multiplier stepdown
 function stepdown!{T<:AbstractFloat}(sortedPValues::Vector{T}, multiplier::Function, n::Integer = length(sortedPValues))
-  stepfun(p::T, i::Int, n::Int) = p*multiplier(i,n)
+  stepfun(p::T, i::Int, n::Int) = p * multiplier(i, n)
   general_stepdown!(sortedPValues, stepfun, n)
   return sortedPValues
 end
 
 function general_stepdown!{T<:AbstractFloat}(sortedPValues::Vector{T}, stepfun::Function, n::Integer = length(sortedPValues))
     sortedPValues[1] = stepfun(sortedPValues[1], 1, n)
     for i in 2:n
-        sortedPValues[i] = max(sortedPValues[i-1], stepfun(sortedPValues[i],i, n))
+        sortedPValues[i] = max(sortedPValues[i-1], stepfun(sortedPValues[i], i, n))
     end
     return sortedPValues
 end
 
+
 function reorder{T<:Number}(values::Vector{T})
     newOrder = sortperm(values)
     oldOrder = sortperm(newOrder)
@@ -47,34 +48,36 @@ function isin{T<:Real}(x::Vector{T}, lower::Real = 0., upper::Real = 1.)
 end
 
 
-function isotonicregression(y::Array{Float64,1},w::Array{Float64,1})
-  #todo: ignore zero weights
-  y=copy(y)
-  w=copy(w)
-  m = length(y)
-  cnts = ones(Int64,m)
-  i = 2
-  # ... not most efficient way but could be fun to (ab)use iterator protocol
-  while (!done(y,i))
-    if y[i]<y[i-1]
-      y[i-1]=(w[i]*y[i]+w[i-1]*y[i-1])/(w[i]+w[i-1])
-      w[i-1]=w[i]+w[i-1]
-      cnts[i-1] += cnts[i]
-      deleteat!(y,i)
-      deleteat!(w,i)
-      deleteat!(cnts,i)
-      i = max(i-2,1)
+function isotonic_regression_reference{T<:AbstractFloat}(y::Vector{T}, w::Vector{T})
+    #todo: ignore zero weights
+    y = copy(y)
+    w = copy(w)
+    m = length(y)
+    cnts = ones(Int64, m)
+    i = 2
+    # ... not most efficient way but could be fun to (ab)use iterator protocol
+    while !done(y, i)
+        if y[i] < y[i-1]
+            y[i-1] = (w[i]*y[i]+w[i-1]*y[i-1])/(w[i]+w[i-1])
+            w[i-1] = w[i]+w[i-1]
+            cnts[i-1] += cnts[i]
+            deleteat!(y, i)
+            deleteat!(w, i)
+            deleteat!(cnts, i)
+            i = max(i-2, 1)
+        end
+        i += 1
     end
-    i += 1
-  end
-  yisotonic = vcat([y[idx]*ones(Float64,cnt) for (idx,cnt) in enumerate(cnts)]...)
+    yisotonic = vcat([y[idx]*ones(Float64, cnt) for (idx, cnt) in enumerate(cnts)]...)
+    return yisotonic
 end
 
-function isotonicregression(y::Array{Float64,1})
-  isotonicregression(y, ones(Float64, length(y)))
+function isotonic_regression_reference{T<:AbstractFloat}(y::Vector{T})
+    isotonic_regression_reference(y, ones(y))
 end
 
-function isotonic_regression(y::Vector{Float64}, weights::Vector{Float64})
+
+function isotonic_regression{T<:AbstractFloat}(y::Vector{T}, weights::Vector{T})
     n = length(y)
     if n <= 1
         return y
@@ -115,12 +118,12 @@ function isotonic_regression(y::Vector{Float64}, weights::Vector{Float64})
     return y
 end
 
-function isotonic_regression(y::Vector{Float64})
+function isotonic_regression{T<:AbstractFloat}(y::Vector{T})
     isotonic_regression(y, ones(y))
 end
 
 
-function grenander(pv::Vector{Float64})
+function grenander{T<:AbstractFloat}(pv::Vector{T})
     pv_sorted = sort(pv)
     ## ecdf that handles duplicated values
     pv_sorted_unique, counts = rle(pv_sorted)
@@ -137,4 +140,3 @@ function grenander(pv::Vector{Float64})
 
     return pv_sorted_unique, f, F
 end
-

test/test-grenander.jl

@@ -38,12 +38,12 @@ using Base.Test
     @testset "Isotonic regression" begin
 
         #pooled adjacent violators example page 10 robertson
-        isotonicregression = MultipleTesting.isotonicregression
-        @test isotonicregression([22.5; 23.333; 20.833; 24.25], [3.0;3.0;3.0;2.0]) ≈ [22.222; 22.222; 22.222; 24.25]
+        isotonic_regression_reference = MultipleTesting.isotonic_regression_reference
+        @test isotonic_regression_reference([22.5; 23.333; 20.833; 24.25], [3.0;3.0;3.0;2.0]) ≈ [22.222; 22.222; 22.222; 24.25]
 
         #if input already ordered, then output should be the same
-        @test isotonicregression([1.0; 2.0; 3.0]) ≈ [1.0; 2.0; 3.0]
-        @test isotonicregression([1., 41., 51., 1., 2., 5., 24.], [1., 2., 3., 4., 5., 6., 7.]) ≈ [1.0, 13.95, 13.95, 13.95, 13.95, 13.95, 24]
+        @test isotonic_regression_reference([1.0; 2.0; 3.0]) ≈ [1.0; 2.0; 3.0]
+        @test isotonic_regression_reference([1., 41., 51., 1., 2., 5., 24.], [1., 2., 3., 4., 5., 6., 7.]) ≈ [1.0, 13.95, 13.95, 13.95, 13.95, 13.95, 24.]
 
         # single value or empty vector remains unchanged
         r = rand(1)