bspec (version 1.6)

bspec: Computing the spectrum's posterior distribution

Description

Derives the posterior distribution of the spectrum of one or several time series, based on data and prior specifications.

Usage

bspec(x, ...)
  # S3 method for default
bspec(x, priorscale=1, priordf=0, intercept=TRUE,
    two.sided=FALSE, ...)

Arguments

x

a time series object of the data to be analysed. May be a univariate (ts object) or multivariate (mts object) time series.

priorscale

either a numeric vector giving the scale parameters of the spectrum's prior distribution; recycled if of length 1.

Or a function of frequency.

priordf

either a numeric vector giving the degrees-of-freedom parameters of the spectrum's prior distribution; recycled if of length 1.

Or a function of frequency.

intercept

a logical flag indicating whether to include the ‘intercept’ (zero frequency) term.

two.sided

a logical flag indicating whether to refer to a one-sided or a two-sided spectrum. In particular affects the interpretation of the prior scale parameters, and sets the default for some methods applied to the resulting bspec object via its two.sided element.

...

currently unused.

Value

A list of class bspec containing the following elements:

freq

a numeric vector giving the (Fourier-) frequencies that the spectral parameters correspond to.

scale

a numeric vector giving the scale parameters of the posterior distributions of the spectral parameters corresponding to the above frequencies. These -internally- always correspond to the one-sided spectrum, regardless of the two.sided flag (see below).

df

a numeric vector giving the degrees-of-freedom parameters of the posterior distributions of the spectral parameters corresponding to the above frequencies.

priorscale

a numeric vector giving the prior scale parameters.

priordf

a numeric vector giving the prior degrees-of-freedom parameters.

datassq

a numeric vector giving the sum-of-squares contributed by the data.

datadf

a numeric vector giving the degrees-of-freedom contributed by the data.

N

the sample size of the original time series.

deltat

the sampling interval of the original time series.

deltaf

the frequency interval of the Fourier-transformed time series.

start

the time of the first observation in the original time series.

call

an object of class call giving the function call that generated the bspec object.

two.sided

a logical flag indicating whether the spectrum is to be interpreted as one-sided or two-sided.

Details

Based on the assumptions of a zero mean and a finite spectrum, the posterior distribution of the (discrete) spectrum is derived. The data are modeled using the Maximum Entropy (Normal) distribution for the above constraints, and based on the prior information about the spectrum specified in terms of the (conjugate) scaled inverse \(\chi^2\) distribution.

For more details, see the references.

References

Roever, C., Meyer, R., Christensen, N. Modelling coloured residual noise in gravitational-wave signal processing. Classical and Quantum Gravity, 28(1):015010, 2011. 10.1088/0264-9381/28/1/015010. See also arXiv preprint 0804.3853.

Roever, C. Degrees-of-freedom estimation in the Student-t noise model. Technical Report LIGO-T1100497, LIGO-Virgo Collaboration, 2011.

Roever, C. A Student-t based filter for robust signal detection. Physical Review D, 84(12):122004, 2011. 10.1103/PhysRevD.84.122004. See also arXiv preprint 1109.0442.

See Also

expectation, quantile.bspec, sample.bspec, ppsample, acf.bspec, spectrum

Examples

# NOT RUN {
# determine spectrum's posterior distribution
# (for noninformative prior):
lhspec <- bspec(lh)
print(lhspec)

# show some more details:
str(lhspec)

# plot 95 percent central intervals and medians:
plot(lhspec)

# draw and plot a sample from posterior distribution:
lines(lhspec$freq, sample(lhspec), type="b", pch=20)

########

# compare the default outputs of "bspec()" and "spectrum()":
bspec1    <- bspec(lh)
spectrum1 <- spectrum(lh, plot=FALSE)
plot(bspec1) 
lines(spectrum1$freq, spectrum1$spec, col="blue")
# (note -among others- the factor 2 difference)

# match the outputs:
# Need to suppress  tapering, padding and de-trending
# (see help for "spec.pgram()"):
spectrum2 <- spectrum(lh, taper=0, fast=FALSE, detrend=FALSE, plot=FALSE)
# Need to drop intercept (zero frequency) term:
bspec2    <- bspec(lh, intercept=FALSE)
# plot the "spectrum()" output:
plot(spectrum2)
# draw the "bspec()" scale parameters, adjusted
# by the corresponding degrees-of-freedom,
# so they correspond to one-sided spectrum:
lines(bspec2$freq, bspec2$scale/bspec2$datadf,
      type="b", col="green", lty="dashed")

########

# handle several time series at once...
data(sunspots)
# extract three 70-year segments:
spots1 <- window(sunspots, 1750, 1819.99)
spots2 <- window(sunspots, 1830, 1899.99)
spots3 <- window(sunspots, 1910, 1979.99)
# align their time scales:
tsp(spots3) <- tsp(spots2) <- tsp(spots1)
# combine to multivariate time series:
spots <- ts.union(spots1, spots2, spots3)
# infer spectrum:
plot(bspec(spots))
# }