Matlab code for estimating linear Markov models

Copyright (C) 2006 Matthew Spencer, School of Biological Sciences, University of Liverpool (m.spencer@liverpool.ac.uk).

Dependencies

Matlab (tested with version 7.3.0.298 (R2006b), The Mathworks Inc, Natick, MA
Matlab optimization toolbox (tested with version 3.1)

--------------------------------------------------------------

To estimate parameters for a continuous-time linear model:

optresults_hom_randinit=gethomlike_randinit(n0,ntrans_cat,tset,fname,nreps,randseed,gethessian)

where:

n0 is a column vector of initial state counts

ntrans_cat is an array of concatenated matrices of transition counts [F1,F2,F3,...]
For a matrix Fk, entry f_{ij} is the number of observations of points in state j at time k-1 that were in state i at time k

tset is the set of times k (excluding zero), that are the ends of each time interval

fname is the name of a file to which the optimization results will be saved

nreps is the number of replicate optimizations to run

randseed (optional) is a seed for the pseudorandom number generator (current time is used if not set)

gethessian (optional) is logical: if true, compute the hessian at the solution (default is not to compute)

The results are written to the structure optresults_hom_randinit, which has fields (for each replicate optimization):

negl is negative log likelihood

p0_est is estimated initial frequencies

Q_est is estimated rate matrix

exitflag indicates exit condition of optimizer fmincon

optoutput is a structure with information about the optimization

grad is the gradient of the likelihood function at the solution

hessian is the value of the Hessian (second derivative) matrix at the
solution, for the transformed parameters used in unconstrained optimization

timetaken is clock time used

init_rates is initial rate estimates

Example:

>> load LVsim_testdata
>> optresults_hom_randinit=gethomlike_randinit(n0,ntrans_cat,tset,fname,nreps,[],gethessian)

fits two replicate optimizations to the data in LVsim_testdata, saves the results to LVsim_testresults.mat, using current clock time as the seed and does not estimate the Hessian.

--------------------------------------------------------------

To estimate parameters for a discrete time-averaged model:

[negl_dtav,Pav,p0_av,nparm_av]=get_dtav_like(n0,ntrans_cat,tset)

n0, ntrans_cat, tset as above

Outputs are:

negl_dtav is negative log likelihood

Pav is discrete time-average transition probability matrix estimate

p0_av is initial frequency estimates

nparm_av is number of parameters

Example:

>> load LVsim_testdata
>> [negl_dtav,Pav,p0_av,nparm_av]=get_dtav_like(n0,ntrans_cat,tset)

--------------------------------------------------------------

To estimate parameters for a discrete saturated model:

[neglobs,p0_saturated,Pset_saturated,nparm_saturated]=LV_saturated_like(n0,ntrans_cat,tset);

n0, ntrans_cat, tset as above

Outputs are:

neglobs is the -log likelihood for the saturated model

p0_saturated is saturated (multinomial) estimates of initial freqs

Pset_saturated is concatenated set of saturated transition probability estimates for each time interval (separate multinomial for each interval)

nparm_saturated is number of parameters in saturated model

Note that each of these functions makes use of other .m files included in the distribution