Multi-layered linear stochastic differential equation inference
This page presents the supplementary material and software used for
multi-layered analysis using vectorial linear
stochastic differential equations. The reason for using stochastic
differential equation is to facilitate studying irregularly spaced
observations like those we find in our study dataset of micro-fossils.
This work was performed in
conjuction with doing the research for and writing the manuscript
"Modelling Evolution by Linear Stochastic Differential Equations with Latent
T. Reitan, T. Schweder and J. Henderiks.
A short powerpoint presentation can be found
here, with a PDF version
here. (Made for the Seventh Conference on
Bayesian Inference in Stochastic Processes in Madrid, Spain,
A simpler presentation of the basic modelling framework and how it
connects to the application can be
found here. A how-to concerning
the program usage, with a short intro to continuous time modelling
is found here.
Originally, this was used for modelling phenotypic
evolution towards a optima that itself could be a stochastic process
depending responding to layers below it (climate or primary optimum).
However, it may well be that the method can be used for other
applications as well, with a little tweaking. The program found here
so far, was geared towards solving the specific problem of Coccolithus
microfossils though, but has lately been made to solve more general
linear continuous time process inference problems.
The program (as well as the previously used underlying program library
Hydrasub) is licenced as
meaning that anyone can modify and use it
for other purposes than the original intent.
The latest version of the analysis program and R package uses only
the Lapack library.
Older versions use the Hydrasub library for
support. This library in turn uses GSL
(Gnu Scientific Library). Plotting programs, such as 'vvgraph' and 'histogramme'
are used for displaying the results. PS: If plotting in the program (rather than the
R package) is to be used, then even the latest version requires the
'vvgraph' and 'histogramme' programs from the Hydrasub library.
However, there are options for not doing plots and for sending plots to file instead.
- supplementary.pdf: A
supplementary to the article "Modelling Evolution by Linear
Stochastic Differential Equations with Latent Layers", which
goes more into detail concerning the algorithmical aspects
of the inference.
- layeranalyzer.cpp. This is the
current source version of the program (replacing "layer_analyzer_complex"
see below). This version does not depend on any library but the
Lapack library (for linear algebra operations), which typically can be
found in an R installation. The program can then typically be
compiled on a Linux machine like this:
"g++ -DMAIN -DENGLISH_LANGUAGE -I/usr/include/R -I/usr/include/R/R_ext -o layeranalyzer layeranalyzer.cpp -lm -llapack". If the Lapack library files or header files
are somewhere else on your computer, you should change the "-I" (which
tells the compiler where to look for header files) and "-L" (which tells
the compiler where to search for library files other than in the
- An R package,
for the same analyzis has been made, which uses the same
underlying C++ code (layeranalyzer.cpp).
It can be installed with the R code:
install.packages("http://folk.uio.no/trondr/R/layeranalyzer_0.1.tar.gz",type="source",verbose=T), and used after invoking library(layeranalyzer).
This has already been tested on multiple platforms. See "help" texts for more
layer_analyzer but is able to deal with feedback loops and complex
eigenvalues (thus making this a tractable way of continuous time
ecological modelling). May run slower on such problems, but should run just
as fast as usual on models without feedback loops. Also able to
deal with non-stationary (exploding) time series, as long as
an initial value treatment is specified. Also, multiple measurements
at the same time possble (as long as they don't belong to the same
site and the same series). Thus phenotypical traits with
allometry can be handled. Replaced the old "layer_analyzer"
as the default analysis program, but has now been replaced by the
program "layeranalyzer" (see above, and sorry for the lack of imagination),
which doesn not use the Hydrasub library and has replaced the use of the
GSL library for the more universally available "lapack" library.
layer_analyzer.H: A generalized
layered analyzer, though has been replaced by an even
more generalized program (now called layer_analyzer_complex),
which allows for feedback loops (something this version does not).
It does Bayesian analysis (and also
a subsequent ML analysis, if that is wanted) on a
layered stochastic differential equation set. More than one
time series can be analyzed at a time, with options for letting
a layer belonging to one series affect a layer belonging to another.
A Perl-script for traversing a
host of models (710 in number)
in order to find the posterior model probabilities and find the most
probable model. Uses 'layer_analyzer'.
This is the aggregated dataset with smoothed standard deviations
analyzed in the paper
"Modelling Evolution by Linear Stochastic Differential Equations
with Latent Layers", collected by Jorijntje Henderiks.
contains the prior specification, used in layered_analyzer.
- default_prior.txt. Defines the
default prior used in the R package.
Various similated test files.
- temp.txt: Three series of water
temperature from the Norwegian Water Resources Administration,
NVE. Some data has been removed
to check the state inference. The full data is found in this file:
temp.txt which can be used for comparison
in layered_analyzer. Prior file:
A bundle of the supplementary text,
the dataset and
the main source codes (layer_analyzer).
- ou_test2.R: R code for comparing
MCMC+importance sampling BML estimation with (semi-)exact
BML (done by discresizing characteristic time and using conjugate
priors for expectancy and variance). Performed on an OU model
with artificial data, which shows that the method can come
really close to the exact BML. The code uses a very simple
MCMC algorithm (RW Metropolis with adaptive burn-in) plus the
importance sampling method itself, both of which can be found in
this file: mcmc.R.
- 23/11-2017: New version of the analysis program simply
called "layeranalyzer" added. This version does not depend on
anything but the Lapack library. All parts of the Hydrasub library
that were needed, has been included in the program code itself.
This should make compilation easier and allow for compilation on
more platforms. R package also added.
- 25/2-2015: Debug on layer_analyzer_complex has been performed
on the option for linear time trends and the option has been tested on
simulated files. The input for multiple files has been shortened
down to avoid using having dual ways of specifying the same options.
Run layer_analyzer_complex without input arguments to see the changes.
- 21/3-2013: The program (layer_analyzer_complex) can now
subtract a periodic
pattern (sine/cosine) from the observations, thus allowing seasonality or
other strict periodic phenomena to be removed. The length of
the period is specified. Multiple periodicities can be specified,
making ti possible to handle arbitrarily complex period patterns.
The coefficients in fron of the sine/cosine terms enter in as parameters.
In addition, some extra debug info was added as options. Also
the parallel tempering was given an extra option, where the
multiplicative steps of the 'temperature' now can be specified.
Also, starting parameters for the MCMC chain can now be set.
A short intro to using the program was made,
here. The into does not show
all the options. Call the program without input in order to see that.
- 21/5-2012: The program layer_analyzer_complex is now
able to deal with multiple correlated series, using an
extra file specifying the correlation in the observations. Thus
allometry in biology is dealt with in the model. Correlations
as between series as well as causal pathways are possible. This
machinery is currently being tested on microfossil data. The
results so far look encouraging.
- 23/3-2012: A large expansion that will later replace
the current layer_analyzer is now made available as a
separate program named layer_analyzer_complex. The program deals
with separately specified feedback loops and also with the resulting
problem of complex eigenvalues.
- 7/2-2012: New code for dealing with pair-wise inter-regional
correlations in a specified layer. (Note that there should be no
more than 6 or at worst 7 sites for the methodology to be
efficient. If there's more, the resources to the initial Monte Carlo
estimation of the probability of drawing positive covariance
structures needs to be adjusted in the code. Either that or a
completely different prior for the correlation structure, like the
inverse-Wishart must be implemented.) Also, an indicator on/off
option for allowing the analysis to switch off inter-regional
correlations for some sites.
- 10/2-2012: Bug fix for 'layer_analyzer' concerning the prior
identification strategies in combination with itnegrated layers.
Also, test code for importance sampling added.
- 17/1-2012: Bug fix on 'layer_analyzer' on
the priors for hierarchical models with
identification restriction. Possibility of switching off
identification restrictions. Realizations from smoothing estimates
as well as simulations based on given parameters also implemented
as options. Possibility of plugging in secondary data sources that
will also be modelled using linear SDEs. The point is to make
'layer_analyzer' a single tool for all practical linear SDE problems.
Up to my home page.
Last modified: Wed Nov 23 12:00:00 CET 2017