gibbs_np: Gibbs sampler for Bayesian nonparametric inference with...

In beyondWhittle: Bayesian Spectral Inference for Time Series

Gibbs sampler for Bayesian nonparametric inference with Whittle likelihood

Description

Obtain samples of the posterior of the Whittle likelihood in conjunction with a Bernstein-Dirichlet prior on the spectral density.

Usage

gibbs_np(
  data,
  Ntotal,
  burnin,
  thin = 1,
  print_interval = 100,
  numerical_thresh = 1e-07,
  M = 1,
  g0.alpha = 1,
  g0.beta = 1,
  k.theta = 0.01,
  tau.alpha = 0.001,
  tau.beta = 0.001,
  kmax = 100 * coars + 500 * (!coars),
  trunc_l = 0.1,
  trunc_r = 0.9,
  coars = F,
  L = max(20, length(data)^(1/3))
)

Arguments

data	numeric vector; NA values are interpreted as missing values and treated as random
Ntotal	total number of iterations to run the Markov chain
burnin	number of initial iterations to be discarded
thin	thinning number (postprocessing)
print_interval	Number of iterations, after which a status is printed to console
numerical_thresh	Lower (numerical pointwise) bound for the spectral density
M	DP base measure constant (> 0)
g0.alpha, g0.beta	parameters of Beta base measure of DP
k.theta	prior parameter for polynomial degree k (propto exp(-k.theta*k*log(k)))
tau.alpha, tau.beta	prior parameters for tau (inverse gamma)
kmax	upper bound for polynomial degree of Bernstein-Dirichlet mixture (can be set to Inf, algorithm is faster with kmax<Inf due to pre-computation of basis functions, but values 500<kmax<Inf are very memory intensive)
trunc_l, trunc_r	left and right truncation of Bernstein polynomial basis functions, 0<=trunc_l<trunc_r<=1
coars	flag indicating whether coarsened or default bernstein polynomials are used (see Appendix E.1 in Ghosal and van der Vaart 2017)
L	truncation parameter of DP in stick breaking representation

Details

Further details can be found in the simulation study section in the references papers.

Value

list containing the following fields:

psd.median, psd.mean	psd estimates: (pointwise) posterior median and mean
psd.p05, psd.p95	pointwise credibility interval
psd.u05, psd.u95	uniform credibility interval
k, tau, V, W	posterior traces of PSD parameters
lpost	trace of log posterior

References

C. Kirch et al. (2018) Beyond Whittle: Nonparametric Correction of a Parametric Likelihood With a Focus on Bayesian Time Series Analysis Bayesian Analysis <doi:10.1214/18-BA1126>

N. Choudhuri et al. (2004) Bayesian Estimation of the Spectral Density of a Time Series JASA <doi:10.1198/016214504000000557>

S. Ghosal and A. van der Vaart (2017) Fundamentals of Nonparametric Bayesian Inference <doi:10.1017/9781139029834>

Examples

## Not run: 

##
## Example 1: Fit the NP model to sunspot data:
##

data <- sqrt(as.numeric(sunspot.year))
data <- data - mean(data)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_np(data=data, Ntotal=10000, burnin=4000, thin=2)

# Plot spectral estimate, credible regions and periodogram on log-scale
plot(mcmc, log=T)


##
## Example 2: Fit the NP model to high-peaked AR(1) data
##

n <- 256
data <- arima.sim(n=n, model=list(ar=0.95)) 
data <- data - mean(data)
omega <- fourier_freq(n)
psd_true <- psd_arma(omega, ar=0.95, ma=numeric(0), sigma2=1)

# If you run the example be aware that this may take several minutes
print("example may take some time to run")
mcmc <- gibbs_np(data=data, Ntotal=10000, burnin=4000, thin=2)

# Compare estimate with true function (green)
plot(mcmc, log=F, pdgrm=F, credib="uniform")
lines(x=omega, y=psd_true, col=3, lwd=2)

# Compute the Integrated Absolute Error (IAE) of posterior median
cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")

## End(Not run)

beyondWhittle documentation built on June 22, 2024, 11:35 a.m.