5 files changed, 129 insertions(+), 45 deletions(-)

M .gitignore
M README.md
M example/basicexample.jl
A example/basicmodel.jl
A example/fitsimulation.jl

@@ 1,5 1,6 @@
 *.csv
 *.png
+*.jld
 !simulationplot.png
 plotsim.gp
 cmocean_balance.pal

@@ 7,14 7,13 @@ nonstationary spatio-temporal modeling of high-frequency monitoring data* (arxiv
 link coming soon). 
 
 See the paper for a discussion of the covariance function itself, and see the
-example file `./example/basicexample.jl` for a quick-start guide. With just a
-few lines of code, you can specify a complicated covariance function
-corresponding to fields like this:
-
-<!---
-![](./simulationplot.png "An example simulation from the
-covariance function created in `./example/basicexample.jl`")
---->
+example file `./example/basicmodel.jl` for an example specification of the
+model, `./example/basicexample.jl` for an example of creating a kernel function
+and then simulating a process with it, and finally `./example/fitsimulation.jl`
+for an example of computing an MLE for the simulated data. With very few lines
+of boilerplate, you can specify  a complicated covariance function corresponding
+to fields that look like this:
+
 <p align="center">
   <img src="https://git.sr.ht/~cgeoga/HalfSpectral.jl/blob/master/simulationplot.png">
 </p>


@@ 23,10 22,10 @@ For the code used to fit models in the paper, see `./examples/paperscripts/`. I
 do not plan to modify these files if this code evolves, so if you want to run
 those scripts, please use the appropriate tagged release.
 
-For estimation, please see my other software for estimation,
-[GPMaxlik.jl](https://git.sr.ht/~cgeoga/GPMaxlik.jl). The examples for setting
-up and performing optimization in that repository are much simpler than the
-examples in `./examples/paperscripts/`.
+For more information about estimation, please see another of my software
+packages, [GPMaxlik.jl](https://git.sr.ht/~cgeoga/GPMaxlik.jl). All the
+necessary code to try simulating and estimating is here, but several more
+complete examples are provided in that repository.
 
 **If you use this software (HalfSpectral.jl), please cite the paper and not
 the software package itself.**

@@ 1,32 1,13 @@
 
 using HalfSpectral, LinearAlgebra
 
-# Marginal spectral density, with the twist of making a few parameters depend on
-# the spatial index.
-function Sf(w,x,p)
-  (scale, rate, smoothness) = p[1:3]
-  rate += log(x)
-  smoothness *= (1+x/10)
-  scale*(rate*sinpi(w)^2 + 1)^(-smoothness-1/2) 
-end
-
-# A coherence function that does happen to be stationary. See the paperscripts
-# model for code implementing the Paciorek-Schervish nonstationary coherence.
-function Cw(w,x,y,p)
-  gamxy = p[4]/(1 + (sinpi(w)/0.3)^2)^4
-  exp(-abs(x-y)/gamxy)
-end
-
-# A simple phase function.
-phase(w,x,y,p) = exp(2.0*pi*im*p[5]*sin(sinpi(w))*(x-y))
-
-# The integrand in the form that HalfSpectral requests.
-integrand(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)*phase(w,x,y,p) 
-
-# Create the kernel for spatial locations 1:40 and 200 unit time measurements.
-const tg     = timegrid(1:200)
-const xpts   = 1:40
-const params = (50.0,100.0,0.5,5.0,0.5) # sample parameter values.
+# Load a basic specification of a nonstationary half-spectral model.
+include("basicmodel.jl")
+
+# Create the kernel for spatial locations 1:30 and 100 unit time measurements.
+const tg     = timegrid(1:100)
+const xpts   = 1:30
+const params = (log(50.0),log(100.0),log(0.5),log(5.0),0.5) # sample parameter values.
 kernel = tftkernel(integrand, tg, xpts, params)
 
 # Build a second time to get a timing example.


@@ 41,17 22,31 @@ println("Building the kernel takes this long")
 # disagrees with the internal one it will rebuild the dictionary of function
 # values. In general, I suggest including the parameter vector for code safety.
 # But to avoid runtime checks, you can use the unchecked kernel(x,y).
-pts = vec(reverse.(collect(Iterators.product(xpts, tg.Xt))))
-K   = Symmetric([kernel(x,y) for x in pts, y in pts])
-sim = reshape(cholesky(K).L*randn(length(pts)), length(xpts), length(tg.Xt))
+pts  = vec(reverse.(collect(Iterators.product(xpts, tg.Xt))))
+K    = Symmetric([kernel(x,y) for x in pts, y in pts])
+simv = cholesky(K).L*randn(length(pts))
+simm = reshape(simv, length(xpts), length(tg.Xt))
+
+# Optionally, save the simulation file.
 
-#= Optionally, save the simulation file:
+#= If you just want a csv file to plot.
 using DelimitedFiles
-writedlm("sim.csv", sim, ',')
+writedlm("sim.csv", simm, ',')
 =#
 
-#= Optionally, visualize it:
+# If you want to save enough output to also fit your simulated data:
+using JLD
+JLD.save("simulation_output.jld",
+         "true_params", params,
+         "pts",   pts,
+         "tgrid", tg,
+         "xpts",  xpts,
+         "data",  simv)
+# =#
+# =#
+
+#= Optionally, visualize it.
 using Plots
-heatmap(sim, xlabel="time", ylabel="altitude")
+heatmap(simm, xlabel="time", ylabel="altitude")
 =#
 

@@ 0,0 1,31 @@
+
+# Marginal spectral density, with the twist of making a few parameters depend on
+# the spatial index.
+function Sf(w,x,p)
+  (scale, rate, smoothness) = (exp(p[1]), exp(p[2]), exp(p[3]))
+  rate += log(x)
+  smoothness *= (1+x/10)
+  scale*(rate*sinpi(w)^2 + 1)^(-smoothness-1/2) 
+end
+
+# A coherence function that does happen to be stationary. See the paperscripts
+# model for code implementing the Paciorek-Schervish nonstationary coherence.
+function Cw(w,x,y,p)
+  gamxy = exp(p[4])/(1 + (sinpi(w)/0.3)^2)^4
+  exp(-abs(x-y)/gamxy)
+end
+
+# A simple phase function and it's derivative with respect to kernel parameters.
+phase(w,x,y,p) = exp(2.0*pi*im*p[5]*sin(sinpi(w))*(x-y))
+dphase(w,x,y,p)= 2.0*pi*im*sin(sinpi(w))*(x-y)*phase(w,x,y,p)
+
+# The integrand in the form that HalfSpectral requests, splitting up the phase
+# part so we can use AD on the real-valued parts, which seems to avoid a lot of
+# finnicky issues with ForwardDiff.jl. This is not strictly necessary, but as
+# you'll see in fitsimulation.jl, writing it in this slightly tricky way means
+# that with just a few more lines of code, you can use very powerful
+# quasi-Newton optimization methods.
+integrand_nophase(w,x,y,p) = sqrt(Sf(w,x,p)*Sf(w,y,p))*Cw(w,x,y,p)
+integrand(w,x,y,p) = integrand_nophase(w,x,y,p)*phase(w,x,y,p) 
+integrand_dphase(w,x,y,p) = integrand_nophase(w,x,y,p)*dphase(w,x,y,p)
+

@@ 0,0 1,58 @@
+
+using HalfSpectral, GPMaxlik, JLD, ForwardDiff
+
+# Load in the model used to simulate the data.
+include("basicmodel.jl")
+
+# A lazy version of kernel derivatives: use AD on the integrand. This would
+# obviously be much faster if you hand-wrote your derivative functions for the
+# integrand, but as you'll see here this is good enough for prototyping.
+function coord_deriv(f,p,j)
+  ForwardDiff.derivative(z->f(vcat(p[1:(j-1)],z,p[(j+1):end])), p[j])
+end
+dfuns = convert(Vector{Function},
+          [(w,x,y,p)->coord_deriv(z->integrand_nophase(w,x,y,z),p,j)*phase(w,y,x,p)
+           for j in 1:4])
+push!(dfuns, integrand_dphase)
+
+# Load the saved data in:
+const pts   = JLD.load("simulation_output.jld", "pts")
+const tgrid = JLD.load("simulation_output.jld", "tgrid")
+const xpts  = JLD.load("simulation_output.jld", "xpts")
+const dat   = JLD.load("simulation_output.jld", "data")
+const tru   = JLD.load("simulation_output.jld", "true_params")
+
+# Choose some initial conditions, then create the kernel function and the
+# derivative kernel functions.
+ini    = vcat(log.([75.0, 75.0, 1.0, 4.0]), 0.75)
+kfun   = tftkernel(integrand, tgrid, xpts, ini)
+dkfuns = [tftkernel(integrand_dj, tgrid, xpts, ini) for integrand_dj in dfuns]  
+
+# Create some SAA vectors for faster derivatives.
+const saa = rand([-1.0, 1.0], length(pts), 100)
+
+# Create objective only and objective+grad+fisher functions for GPMaxlik, using
+# the optimized nll function it provides.
+obj(p) = nll(pts, dat, kfun, dkfuns, p, saa=saa, nll=true, grad=false, fish=false).nll
+function objgradfish(p)
+  nll_call = nll(pts, dat, kfun, dkfuns, p, saa=saa, nll=true, grad=true, fish=true)
+  (nll_call.nll, nll_call.grad, nll_call.fish)
+end
+
+# Perform the optimization using the trust region method (see GPMaxlik.jl for
+# details) using the stochastic expected Fisher matrix as a Hessian
+# approximation. If you want to try BFGS to compare, just change "fish=true" to
+# "fish=false" in line 38 above.
+maxlik_result = trustregion(obj, objgradfish, ini, vrb=true, maxit=100, rtol=1.0e-5)
+
+# Consider the resulting MLE, std approximations from the stochastic expected
+# Fisher matrix, and compare with the true parameters for the simulation.
+mle   = maxlik_result.p
+stds  = map(x->"("*string(x)*")", 
+            round.(sqrt.(diag(inv(maxlik_result.h))).*1.96, digits=2))
+table = hcat(pushfirst!(string.(round.(collect(tru), digits=4)), "true"),
+             pushfirst!(string.(round.(collect(mle), digits=4)), "mle"),
+             pushfirst!(stds, "±95%"))
+display(table)
+
+