Hey all—
Just a quick announcement of a couple little pieces of software that some of you may find useful or interesting.
HalfSpectral.jl is a simple little package that implements half spectral covariance functions in a way that facilitates working with them in the time domain for regularly measured data. I used this code to specify a big 25-parameter nonstationary model that I fitted to some very complex micrometeorological data. With really just a few tens of lines of code, you can specify interesting enough nonstationary models to obtain simulations like this:
GPMaxlik.jl is another simple little package for maximum likelihood estimation for mean-zero Gaussian processes. It provides a relatively optimized function for evaluating the log-likelihood, it’s gradient and expected Fisher information matrix, and nicely stable stochastic estimators for those two quantities that avoid all matrix-matrix arithmetic. There are also some little tricks for using stochastic estimators until you get near the MLE or reach a certain number of iterations, after which the function calls switch to exact methods, so that you can sort of cheaply get closer to your MLE. I also implemented several convenient ways to inject sparsity information about covariance matrices into this function.
The second thing it implements is a trust region-based optimization method that I more or less implemented line-by-line reading Nocedal and Wright (2006). For whatever reason, trust region methods seem less popular in the Julia ecosystem, but in my experience they work much better for maximum likelihood problems than line search methods. I am not an optimizer by any stretch and this code is not sophisticated, but it’s pretty flexible and implemented in such a way that you can easily inject your own print outputs or stopping conditions and so on. This is what I used to fit the huge model described above, and I had a better experience with this code than with IPopt, which is really an excellent piece of software. I interpret this to be evidence that trust region methods really can work noticeably better for maximum likelihood problems.
In both repositories, I implemented pretty complete examples that attempt to demonstrate all of the boilerplate (which is pretty small) required to actually specify a model (for HalfSpectral.jl) and then fit it (for GPMaxlik.jl).