MEM: Spatial model

spm.jl

module spm

using Distributions, LinearAlgebra
import MeasurementErrorModels as MEM

function Pf(s::T, p::T) where T<:AbstractVector
    Dist = (s .-s').^2
    return p[2]*exp.(-Dist./p[1])
end

function gls(Y::T, X::T, Σy::T) where T<:AbstractMatrix
    Xa = [ones(size(X,1)) X]
    B = (Xa'/Σy*Xa) \ (Xa'/Σy*Y)
    Σb = inv(Xa'/Σy*Xa)
    covB = 0.5*(Σb+Σb')
    return MvNormal(B[:], covB)
end

# "Weighted total least squares"
function wtls(Y::AbstractVector, X::T, Ωv::T, ΣA::T) where T<:AbstractMatrix
    n, p  = size(X)
    Xa = [ones(n) X]
    Bml1 = MEM.Bml(Y, X, Ωv, ΣA)
    e = Y - Xa*Bml1
    BI = kron(Bml1[2:end], I(n))
    W =  Ωv + BI'ΣA*BI
    covB = MEM.BBml(e, W, Xa, ΣA, BI)
    return MvNormal(Bml1,  covB)
end



end

An Example:

using LinearAlgebra, Distributions, DataFrames, Random

Random.seed!(123)
n = 200
spc = [1.0:n;] # spatial coordinates
Σεy =  spm.Pf(spc, [10.0, 5.0])
Σεx = spm.Pf(spc, [0.5, 5.0])
τ² = 5.0
εx = rand(MvNormal(zeros(n), Σεx),1)
εy = rand(MvNormal(zeros(n), Σεy),1)
e = rand(Normal(0.0, sqrt(τ²)), n)
# std([εx εy e], dims=1)

b = [1.0;;]
Xstar = rand(Normal(0, sqrt(10.0)), n)
X = Xstar + εx
Y = Xstar*b + εy + e
Σy = Σεy + τ²*I(n) 
Bgls = spm.gls(Y, X, Σy)
Bmem = spm.wtls(Y[:], X, Σy, Σεx)

# The true model parameters (intercept & slope)  are: (0.0, 1.0)
B = mean.([Bgls, Bmem])
Db = DataFrame([(b0=b[1],b1=b[2]) for b in B])
Db.method = ["GLS", "MEM"]
Db

Row	b0	b1	method
1	0.256495	0.63151	GLS
2	0.326631	0.981606	MEM