bspec {bspec} R Documentation

## 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, ...)
## Default S3 method:
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.

### 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-squared distribution.

For more details, see the references.

### 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.

### Author(s)

Christian Roever, christian.roever@med.uni-goettingen.de

### References

Roever, C., Meyer, R., Christensen, N. Modelling coloured residual noise in gravitational-wave signal processing. Classical and Quantum Gravity, 28(1):015010, 2011. 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. See also arXiv preprint 1109.0442.

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

### Examples

```# 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:
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))
```

[Package bspec version 1.5 Index]