BaSAR.nest: BSA with nested sampling
In JIC-CSB/BaSAR: Bayesian Spectrum Analysis in R

Description Usage Arguments Details Value Author(s) References Examples

View source: R/BaSAR.R

BSA using nested sampling to calculate the evidence.

1	BaSAR.nest(data, start, stop, nsamples, nbackg, tpoints, nposts)

`data`	data as a 1-dimensional vector
`start`	lower limit of period of interest, in seconds
`stop`	upper limit of period of interest, in seconds
`nsamples`	number of samples within the interval start-stop
`nbackg`	number of background functions to be added to the model
`tpoints`	vector of time points, in seconds
`nposts`	number of samples to draw from the posterior

Nested sampling is a variant of Markov chain Monte Carlo (MCMC), and repeatedly samples the posterior by ranking the sample points by their likelihoods. This function will return the evidence, which can be used to directly evaluate models. A model with a higher evidence should be preferred.

BaSAR.nest should be preferred over the functions BaSAR.auto and BaSAR.modelratio when there are multiple frequencies present in the data. This is because the model comparison in the mentioned functions is based on a local integration around the point of maximum likelihood, and this approximation is not valid when there are several high peaks in the posterior. Calculating the evidences and comparing those is then a better approach.

A list containing:

`samples`	1D vector of omega
`weights`	1D vector of posterior samples, corresponding to omega
`logZ`	the evidence of the sampled posterior
`logZerror`	error of logZ
`momega`	mean omega estimate
`stomega`	standard deviation of mean omega

Emma Granqvist, Matthew Hartley and Richard J Morris.

Sivia, D. S. and Skilling, J. (2006) Data analysis: a Bayesian tutorial. 2nd Edition. Oxford: Oxford science publications. Oxford University Press.

Granqvist, E., Oldroyd, G. E. and Morris, R. J. (2011) Automated Bayesian model development for frequency detection in biological time series. BMC Syst Biol 5, 97.
http://dx.doi.org/10.1186/1752-0509-5-97

# Create time series with omega = 0.5 and 0.3, plus noise
tpoints = seq(from=1, to=200, length=200)
dpoints = sin(0.5 * tpoints) +sin(0.3 * tpoints) + 0.1 * rnorm(200, 0, 1)
# Plot time series
plot(dpoints, type="l", col="blue", xlab="t", ylab="d(t)")
# Run BaSAR with nested sampling to estimate mean omega and evidence
r = BaSAR.nest(dpoints, 6, 600, 100, 0, tpoints, 500)
# Plot samples drawn from probability distribution
plot(r$samples,exp(r$weight),type="h",pch=20, col="red", ylab="",
xlab=expression(omega),xlim=c(0,1),ylim=c(0,max(exp(r$weight))))
# Note that two peaks appear, omega 0.3 and 0.5