# 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="")
# }
Run the code above in your browser using DataCamp Workspace