Implement LWR

src/idw.jl

@@ -72,9 +72,12 @@ function weights(fitted::FittedIDW, uₒ)
   e = fitted.model.exponent
   δ = fitted.model.distance
   d = fitted.state.data
+  Ω = domain(d)
 
   xₒ = coordinates(uₒ)
-  x(i) = coordinates(centroid(domain(d), i))
+  x(i) = coordinates(centroid(Ω, i))
 
-  [1 / δ(xₒ, x(i)) ^ e for i in 1:nrow(d)]
+  λ(i) = 1 / δ(xₒ, x(i)) ^ e
+
+  map(λ, 1:nelements(Ω))
 end

src/lwr.jl

@@ -3,7 +3,7 @@
 # ------------------------------------------------------------------
 
 """
-    LWR(distance=Euclidean(), weightfun=h -> exp(-3 * h^2))
+    LWR(weightfun=h -> exp(-3 * h^2), distance=Euclidean())
 
 The locally weighted regression (a.k.a. LOESS) model introduced by
 Cleveland 1979. It is the most natural generalization of [`IDW`](@ref)
@@ -18,8 +18,76 @@ distance-based weights.
 * Cleveland & Grosse 1991. [Computational methods for local
   regression](https://link.springer.com/article/10.1007/BF01890836)
 """
-struct LWR{D,F} <: GeoStatsModel
-  distance::D
+struct LWR{F,D} <: GeoStatsModel
   weightfun::F
+  distance::D
+end
+
+LWR(weightfun) = LWR(weightfun, Euclidean())
+LWR() = LWR(h -> exp(-3 * h ^ 2))
+
+struct LWRState{D<:AbstractGeoTable}
+  data::D
+end
+
+struct FittedLWR{M<:LWR,S<:LWRState}
+  model::M
+  state::S
+end
+
+status(fitted::FittedLWR) = true
+
+#--------------
+# FITTING STEP
+#--------------
+
+function fit(model::LWR, data)
+  # record state
+  state = LWRState(data)
+
+  # return fitted model
+  FittedLWR(model, state)
 end
 
+#-----------------
+# PREDICTION STEP
+#-----------------
+
+predict(fitted::FittedLWR, var, uₒ) = first(lwr(fitted, var, uₒ))
+
+function predictprob(fitted::FittedLWR, var, uₒ)
+  μ, σ² = lwr(fitted, var, uₒ)
+  Normal(μ, √σ²)
+end
+
+function lwr(fitted::FittedLWR, var, uₒ)
+  w = fitted.model.weightfun
+  δ = fitted.model.distance
+  d = fitted.state.data
+  c = Tables.columns(values(d))
+  z = Tables.getcolumn(c, var)
+
+  Ω = domain(d)
+  n = nelements(Ω)
+
+  x(i) = coordinates(centroid(Ω, i))
+
+  # fit step
+  X = mapreduce(x, hcat, 1:n)
+  Xₗ = [ones(eltype(X), n) X']
+  zₗ = z
+
+  # predict step
+  xₒ = coordinates(uₒ)
+  δs = map(i -> δ(xₒ, x(i)), 1:n)
+  δs .= δs ./ maximum(δs)
+  Wₗ = Diagonal(w.(δs))
+  θₗ = Xₗ' * Wₗ * Xₗ \ Xₗ' * Wₗ * zₗ
+  xₗ = [one(eltype(xₒ)); xₒ]
+  rₗ = Wₗ * Xₗ * (Xₗ' * Wₗ * Xₗ \ xₗ)
+
+  μ = θₗ ⋅ xₗ
+  σ² = norm(rₗ)
+
+  μ, σ²
+end

test/lwr.jl

@@ -0,0 +1,2 @@
+@testset "LWR" begin
+end

test/runtests.jl

@@ -9,7 +9,7 @@ using Statistics
 using Test, Random
 
 # list of tests
-testfiles = ["krig.jl", "idw.jl"]
+testfiles = ["krig.jl", "idw.jl", "lwr.jl"]
 
 @testset "GeoStatsModels.jl" begin
   for testfile in testfiles