add approximate LAD

Author: jbytecode

CHANGELOG.md

@@ -1,4 +1,5 @@
 # Upcoming Release 
+- Add exact argument for LAD. If exact is true then the linear programming based exact solution is found. Otherwise, a GA based search is performed to yield approximate solutions. 
 
 # v0.8.19
 - Update Satman(2013) algorithm

docs/make.jl

@@ -2,25 +2,20 @@ using Documenter, LinRegOutliers
 
 
 makedocs(
-         format = Documenter.HTML(
-                                  prettyurls = get(ENV, "CI", nothing) == "true",
-                                  collapselevel = 2,
-                                  # assets = ["assets/favicon.ico", "assets/extra_styles.css"],
-                                 ),
-         sitename="LinRegOutliers",
-         authors = "Mehmet Hakan Satman <mhsatman@gmail.com>, Shreesh Adiga <16567adigashreesh@gmail.com>, Guillermo Angeris <angeris@stanford.edu>, Emre Akadal <emre.akadal@istanbul.edu.tr>",
-         pages = [
-                  "Algorithms" => "algorithms.md",
-                  "Diagnostics" =>  "diagnostics.md",
-                  "Types" => "types.md",
-                  "Datasets" => "datasets.md"
-                 ]
-        )
-
-
-deploydocs(
-   repo = "github.com/jbytecode/LinRegOutliers",
+    format = Documenter.HTML(
+        prettyurls = get(ENV, "CI", nothing) == "true",
+        collapselevel = 2,
+        # assets = ["assets/favicon.ico", "assets/extra_styles.css"],
+    ),
+    sitename = "LinRegOutliers",
+    authors = "Mehmet Hakan Satman <mhsatman@gmail.com>, Shreesh Adiga <16567adigashreesh@gmail.com>, Guillermo Angeris <angeris@stanford.edu>, Emre Akadal <emre.akadal@istanbul.edu.tr>",
+    pages = [
+        "Algorithms" => "algorithms.md",
+        "Diagnostics" => "diagnostics.md",
+        "Types" => "types.md",
+        "Datasets" => "datasets.md",
+    ],
 )
 
 
-
+deploydocs(repo = "github.com/jbytecode/LinRegOutliers")

src/LinRegOutliers.jl

@@ -20,6 +20,19 @@ export find_minimum_nonzero
 export @extractRegressionSetting
 
 
+# Compact genetic algorithm
+include("cga.jl")
+import .CGA: cga
+
+# Hooke-Jeeves algorithm
+include("hj.jl")
+# The function hj() is not exported.
+
+# Genetic Algorithm
+include("ga.jl")
+import .GA: ga, RealChromosome
+
+
 # Predefined datasets used in outlier detection literature
 include("data.jl")
 import .DataSets: phones, hbk, stackloss
@@ -129,14 +142,6 @@ import .Hadi94: hadi1994
 include("dataimage.jl")
 import .DataImage: dataimage
 
-# Compact genetic algorithm
-include("cga.jl")
-import .CGA: cga
-
-
-# Genetic Algorithm
-include("ga.jl")
-import .GA: ga, RealChromosome
 
 # Modified and original Satman (2012) algorithms
 include("gwlts.jl")

src/asm2000.jl

@@ -7,7 +7,8 @@ import LinearAlgebra: det
 import Clustering: Hclust, hclust, cutree
 
 
-import ..Basis: RegressionSetting, @extractRegressionSetting, designMatrix, responseVector, zstandardize
+import ..Basis:
+    RegressionSetting, @extractRegressionSetting, designMatrix, responseVector, zstandardize
 import ..Diagnostics: mahalanobisSquaredBetweenPairs
 import ..LTS: lts
 import ..SMR98: majona
@@ -65,7 +66,7 @@ function asm2000(X::Array{Float64,2}, y::Array{Float64,1})::Dict
     #stdfit = standardize(ZScoreTransform, predicteds, dims = 1)
     stdres = zstandardize(resids)
     stdfit = zstandardize(predicteds)
-    
+
     pairs = hcat(stdfit, stdres)
 
     pairs = hcat(resids, predicteds)

src/dataimage.jl

@@ -42,7 +42,10 @@ julia> Plots.plot(di)
 Marchette, David J., and Jeffrey L. Solka. "Using data images for outlier detection." 
 Computational Statistics & Data Analysis 43.4 (2003): 541-552.
 """
-function dataimage(dataMatrix::Array{Float64,2}; distance = :mahalanobis)::Array{RGB{Float64}, 2}
+function dataimage(
+    dataMatrix::Array{Float64,2};
+    distance = :mahalanobis,
+)::Array{RGB{Float64},2}
     d = nothing
     if distance == :mahalanobis
         d = mahalanobisSquaredBetweenPairs(dataMatrix)
@@ -55,7 +58,7 @@ function dataimage(dataMatrix::Array{Float64,2}; distance = :mahalanobis)::Array
     end
     colours = 1.0 .- d / maximum(d)
     n, _ = size(d)
-    colormatrix = Array{RGB{Float64}, 2}(undef, n, n)
+    colormatrix = Array{RGB{Float64},2}(undef, n, n)
     for i = 1:n
         for j = 1:n
             @inbounds colormatrix[i, j] = RGB(colours[i, j])

src/diagnostics.jl

@@ -55,7 +55,7 @@ end
 function mahalanobisSquaredBetweenPairs(pairs::Matrix; covmatrix = nothing)
     n, _ = size(pairs)
     newmat = zeros(Float64, n, n)
-    if isnothing(covmatrix) 
+    if isnothing(covmatrix)
         covmatrix = cov(pairs)
     end
     try
@@ -194,7 +194,7 @@ end
 
 function dffit(X::Array{Float64,2}, y::Array{Float64,1}, i::Int)::Float64
     n, _ = size(X)
-    indices = [j for j in 1:n if i != j]
+    indices = [j for j = 1:n if i != j]
     olsfull = ols(X, y)
     Xsub = X[indices, :]
     ysub = y[indices]
@@ -423,7 +423,7 @@ end
 
 function jacknifedS(X::Array{Float64,2}, y::Array{Float64,1}, k::Int)::Float64
     n, p = size(X)
-    indices = [i for i in 1:n if i != k]
+    indices = [i for i = 1:n if i != k]
     Xsub = X[indices, :]
     ysub = y[indices]
     olsreg = ols(Xsub, ysub)

src/hj.jl

@@ -0,0 +1,49 @@
+module HookeJeeves
+
+export hj
+
+function mutate(par, p, d)
+    newpar = copy(par)
+    newpar[p] += d
+    return newpar
+end
+
+function hj(
+    f::FType,
+    par::Vector{Float64};
+    maxiter = 1000,
+    startstep = 5.0,
+    endstep = 0.0001,
+) where {FType<:Function}
+    p = length(par)
+    currentstep = startstep
+    iter::Int64 = 0
+    while iter < maxiter
+        fold = f(par)
+        fnow = fold
+        for currentp = 1:p
+            mutateleft = mutate(par, currentp, -currentstep)
+            fleft = f(mutateleft)
+            mutateright = mutate(par, currentp, currentstep)
+            fright = f(mutateright)
+            if fleft < fold
+                par = mutateleft
+                fnow = fleft
+            elseif fright < fold
+                par = mutateright
+                fnow = fright
+            end
+        end
+        if fold <= fnow
+            currentstep /= 2
+        end
+        if currentstep < endstep
+            break
+        end
+        iter += 1
+    end
+
+    return Dict("par" => par, "iter" => iter, "step" => currentstep)
+end
+
+end # end of module

src/lad.jl

@@ -8,14 +8,18 @@ using GLPK
 import ..Basis:
     RegressionSetting, @extractRegressionSetting, designMatrix, responseVector, applyColumns
 
+import ..HookeJeeves: hj
+import ..GA: ga
+
 """
 
-    lad(setting)
+    lad(setting; exact = true)
 
 Perform Least Absolute Deviations regression for a given regression setting.
 
 # Arguments
 - `setting::RegressionSetting`: RegressionSetting object with a formula and dataset.
+- `exact::Bool`: If true, use exact LAD regression. If false, estimate LAD regression parameters using GA. Default is true.
 
 # Description 
 The LAD estimator searches for regression the parameters estimates that minimize the sum of absolute residuals.
@@ -51,23 +55,52 @@ Dict{Any,Any} with 2 entries:
 ```
 
 """
-function lad(setting::RegressionSetting)
+function lad(setting::RegressionSetting; exact::Bool = true)
     X, y = @extractRegressionSetting setting
-    return lad(X, y)
+    return lad(X, y, exact = exact)
 end
 
 
 """
 
-    lad(X, y)
+    lad(X, y, exact = true)
 
 Perform Least Absolute Deviations regression for a given regression setting.
 
 # Arguments
 - `X::Array{Float64, 2}`: Design matrix of the linear model.
 - `y::Array{Float64, 1}`: Response vector of the linear model.
+- `exact::Bool`: If true, use exact LAD regression. If false, estimate LAD regression parameters using GA. Default is true.
 """
-function lad(X::Array{Float64,2}, y::Array{Float64,1})
+function lad(X::Array{Float64,2}, y::Array{Float64,1}; exact::Bool = true)
+    if exact
+        return lad_exact(X, y)
+    else
+        return lad_approx(X, y)
+    end
+end
+
+function lad_approx(X::Array{Float64,2}, y::Array{Float64,1})
+    n, p = size(X)
+
+    mins = ones(Float64, p) * 10^6 * (-1.0)
+    maxs = ones(Float64, p) * 10^6
+    popsize = 100
+
+    function fcost(par)
+        return sum(abs.(y .- X * par))
+    end
+
+    garesult = ga(popsize, p, fcost, mins, maxs, 0.90, 0.05, 1, p * 1000)
+    best = garesult[1]
+    hookejeevesresult =
+        hj(fcost, best.genes, maxiter = 109000, startstep = 10.0, endstep = 0.000001)
+    betas = hookejeevesresult["par"]
+    result = Dict("betas" => betas, "residuals" => y .- X * betas)
+    return result
+end
+
+function lad_exact(X::Array{Float64,2}, y::Array{Float64,1})
     n, p = size(X)
 
     m = JuMP.Model(GLPK.Optimizer)

src/py95.jl

@@ -148,7 +148,7 @@ function jacknifedS(
     omittedIndices::Array{Int,1},
 )::Float64
     n, p = size(X)
-    indices = [i for i in 1:n if !(i in omittedIndices)]
+    indices = [i for i = 1:n if !(i in omittedIndices)]
     olsreg = ols(X[indices, :], y[indices])
     e = residuals(olsreg)
     s = sqrt(sum(e .^ 2.0) / (n - p - length(omittedIndices)))

src/satman2013.jl

@@ -116,7 +116,7 @@ function satman2013(X::Array{Float64,2}, y::Array{Float64,1})
 
     result = Dict()
     result["outliers"] = outlierset
-    result["betas"] = betas 
+    result["betas"] = betas
     result["residuals"] = resids
 
     return result