Example

using SpatialBernoulli, Optimization,OptimizationOptimJL, Optim
using LineSearches
using Random
import ForwardDiff
rng = MersenneTwister(1234)

Random.MersenneTwister(1234)

Parameters

# define locations in the unit square
my_locations = vcat(([x y] for x in 0:0.2:1 for y in 0:0.2:1)...)
nlocs = length(my_locations[:, 1])


my_distance = [sqrt(sum(abs2, my_locations[i, :] - my_locations[j, :])) for i in axes(my_locations, 1), j in axes(my_locations, 1)]

# randomly generate a SpatialBernoulli
my_λ = rand(nlocs)# [0.5+ 0.1*i for i in 1:nlocs]
my_range = 0.3
my_sill = 1.0
my_order = 1 / 2

0.5

Create the distribution

d = SB(my_range, my_sill, my_order, my_λ, my_distance)

SB{Float64, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}}(
range: 0.3
sill: 1.0
order: 0.5
λ: [0.7174994867957878, 0.3972417661909725, 0.3461510449044749, 0.37993161041802515, 0.5986453330644911, 0.10722690456674966, 0.25055341227188477, 0.1477507709110385, 0.7036252503113573, 0.895325334352325  …  0.4103701221169189, 0.5292063478079192, 0.6449099381947051, 0.050251835552438284, 0.8812432059960263, 0.25778782119817445, 0.6654845820454082, 0.9928346220494833, 0.15601188066176108, 0.49455759859301895]
h: [0.0 0.2 … 1.2806248474865698 1.4142135623730951; 0.2 0.0 … 1.1661903789690602 1.2806248474865698; … ; 1.2806248474865698 1.1661903789690602 … 0.0 0.19999999999999996; 1.4142135623730951 1.2806248474865698 … 0.19999999999999996 0.0]
ΣU: [1.0 0.513417119032592 … 0.013999278427025617 0.008968424961456372; 0.513417119032592 1.0 … 0.020500597165652824 0.013999278427025617; … ; 0.013999278427025617 0.020500597165652824 … 1.0 0.5134171190325921; 0.008968424961456372 0.013999278427025617 … 0.5134171190325921 1.0]
)

use it like any Distributions.jl distribution

Random draws

n = 2000
y=rand(rng,d,n)

using CairoMakie
yplot=y[:,2]
begin
fig = Figure()
ax = Axis(fig[1, 1])

scatter!(
    ax,
    my_locations[yplot .== 0, 1],
    my_locations[yplot .== 0, 2];
    color = :black,
    label = "y = 0"
)

scatter!(
    ax,
    my_locations[yplot .== 1, 1],
    my_locations[yplot .== 1, 2];
    color = :white,
    strokecolor = :black,
    label = "y = 1"
)

axislegend(ax)
fig
end

Fitting

init_range = 0.5
init_order = 0.5
init_lambda = fill(0.4, nlocs)
init_d = SB(init_range, 1.0, init_order, init_lambda, my_distance)

tdist = maximum(my_distance) / 1
wp = 1.0 .* (my_distance .< tdist)

@timed sol3 = fit_mle_vfast(init_d, y, wp; solver = Optim.LBFGS(
    linesearch = LineSearches.BackTracking()
), order=my_order, return_sol=true, maxiters = 2000)

(value = (SB{Float64, Float64, Float64, Vector{Float64}, Matrix{Float64}, Matrix{Float64}}(
range: 0.3133385271791633
sill: 1.0
order: 0.5
λ: [0.706, 0.3835, 0.354, 0.369, 0.5925, 0.11, 0.2555, 0.146, 0.716, 0.8885  …  0.4245, 0.537, 0.6445, 0.051, 0.887, 0.2725, 0.677, 0.9885, 0.144, 0.4835]
h: [0.0 0.2 … 1.2806248474865698 1.4142135623730951; 0.2 0.0 … 1.1661903789690602 1.2806248474865698; … ; 1.2806248474865698 1.1661903789690602 … 0.0 0.19999999999999996; 1.4142135623730951 1.2806248474865698 … 0.19999999999999996 0.0]
ΣU: [1.0 0.5281962957965838 … 0.01678897546999861 0.010961427875308426; 0.5281962957965838 1.0 … 0.024189841785259934 0.01678897546999861; … ; 0.01678897546999861 0.024189841785259934 … 1.0 0.5281962957965839; 0.010961427875308426 0.01678897546999861 … 0.5281962957965839 1.0]
)
, retcode: Success
u: [-1.1604711164198667]
Final objective value:     2.619638248827536e6
), time = 5.269610731, bytes = 516681032, gctime = 0.06159799, gcstats = Base.GC_Diff(516681032, 8427, 0, 10254536, 66, 6125, 61597990, 1, 0), lock_conflicts = 0, compile_time = 4.552297307, recompile_time = 0.049670849)