1 Million data points

zigzagregression.jl

using ZigZagBoomerang
using Distributions
using SparseArrays
using LinearAlgebra
using Random
Random.seed!(1)
# Million data points, 5 parameters + intercept
N = 1000000
k = 5;
X = randn(N,k);

β = randn(k);

y = X * β + randn(N);


# unbiased estimate of partial derivative of negative log-likelihood with respect to β[k]
# with control variate
function ∇ϕkhat(β, k, samples, X, y, μ) 
    s = 0.0 # ∇priork(β, k) - ∇ϕk(μ, k) 
    for _ in 1:samples
        i = rand(1:length(y))
        s += length(y)/samples*(-X[i,k]'*(y[i]  - dot(X[i,:], β))) 
        s -= length(y)/samples*(-X[i,k]'*(y[i]  - dot(X[i,:], μ)))
    end
    s
end

s = 20 # number of observations per batch

t0 = 0.0
μ = X\y # Perfect control variate: ∇ϕ(μ) = 0
x0 = μ
θ0 = rand([-1.0,1.0], k)
Γ = (sparse((X'*X)))
c = 1*ones(k)
Z = ZigZag(Γ, μ)
T = 1.0
dt = 0.0001
 
@time trace, (tT, xT, θT), (acc, num), re = pdmp(∇ϕkhat, t0, x0, θ0, T, c, Z, s, X, y, μ, adapt=true)

tsd, xsd = ZigZagBoomerang.sep(collect(discretize(trace, dt))) # equidistant for averages
ts, xs = ZigZagBoomerang.sep(collect(trace)) # not equistance for plotting

@show acc, num, acc/num
display(round.(c))
#display(Matrix([Γ inv(cov(xs))]))
display([β mean(xsd)])

using Makie
lines(getindex.(xs,1), getindex.(xs,2))