version 0.18.1

DESCRIPTION

@@ -1,8 +1,8 @@
 Package: beyondWhittle
 Type: Package
 Title: Bayesian Spectral Inference for Stationary Time Series
-Version: 0.17
-Date: 2017-01-19
+Version: 0.18.1
+Date: 2017-04-07
 Authors@R: c(
     person("Alexander", "Meier", email="alexander.meier@ovgu.de", role=c("aut", "cre")),
     person("Claudia", "Kirch", email="claudia.kirch@ovgu.de", role="aut"),
@@ -11,14 +11,14 @@ Authors@R: c(
 Maintainer: Alexander Meier <alexander.meier@ovgu.de>
 Description: Implementations of a Bayesian parametric (autoregressive), a Bayesian nonparametric (Whittle likelihood with Bernstein-Dirichlet prior) and a Bayesian semiparametric (autoregressive likelihood with Bernstein-Dirichlet correction) procedure are provided. The work is based on the corrected parametric likelihood by C. Kirch et al (2017) <arXiv:1701.04846>. It was supported by DFG grant KI 1443/3-1.
 License: GPL (>= 3)
-Imports: ltsa (>= 1.4.6), Rcpp (>= 0.12.7)
+Imports: ltsa (>= 1.4.6), Rcpp (>= 0.12.5)
 LinkingTo: Rcpp
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1
 NeedsCompilation: yes
-Packaged: 2017-01-23 01:48:29 UTC; alexander
+Packaged: 2017-04-07 06:08:55 UTC; alexander
 Author: Alexander Meier [aut, cre],
   Claudia Kirch [aut],
   Matthew C. Edwards [aut],
   Renate Meyer [aut]
 Repository: CRAN
-Date/Publication: 2017-01-23 08:21:26
+Date/Publication: 2017-04-07 08:01:18 UTC

MD5

@@ -1,64 +1,65 @@
-f3389f6516ca797a11103e8c1c7505e7 *DESCRIPTION
-27b30854a0cbf00bd8eea37bc45c89f1 *NAMESPACE
+531c3064a617cfbbadaae6927eb1b864 *DESCRIPTION
+d82819f85b97c7997b9ef7ec0ea7faa0 *NAMESPACE
 39981f1f4f6a547071cf8c29c5ee81d7 *R/RcppExports.R
-7a81c4d46b7ac78698d4bed8948d2105 *R/gibbs_AR.R
-53a43e07179cf3829932e8d646406d4e *R/gibbs_NP.R
-bc2ebb06c513d249697a5dbb372e7140 *R/gibbs_NPC.R
-fa02bc1b5820b81042df00cc869e6563 *R/intern__gibbs_core.R
+8477a9b8ac6b661a0eea312fb208cc1b *R/gibbs_AR.R
+306ad05c319814c996ef02f2fa8cb12f *R/gibbs_NP.R
+aafbbeba51fbd9b6fdc405207d80d4bd *R/gibbs_NPC.R
+b2c266f712f3e5b05f9e58fd56739ca0 *R/intern__gibbs_core.R
 20074a769f4cd4b025e416da990f851b *R/intern__gibbs_pacf.R
-3fce7258a7ef79be3e36fe7edc838450 *R/intern__gibbs_toggle_ar.R
-28d07cdf25ef4ac99a1cf751a499b192 *R/intern__gibbs_util.R
-3526f99df76a510b27ee8e47100d151c *man/Cn.Rd
-0051fcd188f28b47bf291c72f2b56727 *man/acceptanceRate.Rd
-7adfbe557a106ea28f91e1297c8c7290 *man/acvMatrix.Rd
-bacbd124c519a67ca8ed19cd4c2f6bb9 *man/ar_lik.Rd
-c5dfd8180ab15c3ed62b20c284641edf *man/ar_screeType.Rd
-2541c3c8236f74d62d23ec19528dba17 *man/arma_conditional.Rd
-5b06a0e7112ab018438a850ea5d2701b *man/arma_partial.Rd
+736a8ea99e745f73563ab190ab155a6e *R/intern__gibbs_toggle_ar.R
+9fb92549a0e95cacbf63ddc62bd2bff2 *R/intern__gibbs_util.R
+780994827bbfdccb13cadc9e2ad5e5e0 *man/Cn.Rd
+fba22e29ae622ed20351c68b3ab8512a *man/acceptanceRate.Rd
+c917c20e67c25dd899dbf8ce77d7ee8c *man/acvMatrix.Rd
+bf32b0fee7d44800d1624ce205382be9 *man/ar_lik.Rd
+2474fef14e00389b22490104a9410080 *man/ar_screeType.Rd
+5c702345f3bbb8dd5b92ad3ac3921cfb *man/arma_conditional.Rd
+b32fd363e0d747d8e0de0ad1c87b743c *man/arma_partial.Rd
 13d4af3637b07101e1f21ac66546c19c *man/beyondWhittle-package.Rd
-e0721587c62f9d56268d4ea446fb4b1d *man/densityMixture.Rd
-1946ffc6729619f87bb679bad15133ed *man/dtex.kurt.Rd
-290f0ffa83c97fa1b2a3724cd0b0cba0 *man/fast_ft.Rd
-fdc1e1fab7e8d1400837230fedcab6de *man/fast_ift.Rd
-847a3d0d958adef11b0a0247454caa13 *man/fast_mean.Rd
-dde335dce7870827f885139aa6f0203e *man/filenameMCMC.Rd
-954d0ae714b92d36ca66ff7d2eaff932 *man/genEpsARC.Rd
-dada049e9a35c73002d2222ff31d8fb3 *man/genEpsARMAC.Rd
-66497db6383eac3d2c4e770d445a7d26 *man/genEpsMAC.Rd
-e6ef91cfb81a4374d9e28003fc58c57a *man/generalizedGaussian.alpha.Rd
-7898f8892696c4b5ddf1a6a474f7e01e *man/generalizedGaussian.kurtosis.Rd
-7631cb3023dda92637d87b530cbfd8db *man/gibbs_AR.Rd
-2bc6e6c77593be81086290dcf8589e32 *man/gibbs_NP.Rd
-03e53d59ea853403d3cdb8fe39d65940 *man/gibbs_NPC.Rd
-e6d29aabc98d511cb300cf50531a1fff *man/gibbs_pacf.Rd
-d050215a0f67b7cb141681453577516a *man/gibbs_toggle_ar.Rd
-75f82121e9337874195a25c508398a45 *man/l_generalizedGaussian.Rd
-d20c4b78b66e58ee476996e6b0f48ad3 *man/llike.Rd
-ad8037b8f0330cd7e9c24efe28a01445 *man/llike_full.Rd
-d0937fffd0d5fea85f8e60925d46aaa8 *man/llike_pacf.Rd
-45041e8e2d1750aa18d62dc03181feb0 *man/logfuller.Rd
-ae3d363fc4c76effb3c560752c6f8506 *man/logrosenthal.Rd
-6f582dcb446659e738b9905bfd217550 *man/lpost.Rd
-6cf983cade9ff24cdb4f44a49fd6e439 *man/lpost_pacf.Rd
-af0eb6882a74148674eede5e91c53924 *man/lprior.Rd
-09119d27afaa7d8dbcc31eb133126337 *man/lprior_pacf.Rd
-d7711178f75248e24861aa5c28b87686 *man/mixtureWeight.Rd
-69314480f4771ec631d5c8c9679012e8 *man/nll_generalizedGaussian.Rd
-296ddff67e8297e41a585b5c12355941 *man/nll_norm.Rd
-04718bfa5e30c2db45556af066bf1e41 *man/nll_norm_unnormalized.Rd
-ebfe73969242d755a93693598824d452 *man/nll_t.Rd
-162c30142a2c7f1588cda7ed14cf81fb *man/omegaFreq.Rd
-03e33272ae36ac5ddb81595e55b2a591 *man/pFromV.Rd
-0d5a112857d269555bdefb566aec5520 *man/pacf2AR.Rd
-1b284aca0e8a922ffdc461a0b838f6c1 *man/pacf2ARacv.Rd
-f0f704b31bc5f49a98a8ad18b9b9d81b *man/pacfToAR.Rd
-55861392de22ca7ad784601a5365b6f4 *man/plotMCMC.Rd
-d3f0bb9046be2dedd9c575af50d4d365 *man/plotPsdEstimate.Rd
-a27e0efbfa19fe879dce7b68ebcf39bd *man/psd_arma.Rd
-378a552e5fdf90ada5462c871c0843ed *man/reduceMemoryStorageMCMC.Rd
-0d7366d7d102d3c725f46d6cec4408e9 *man/se_kurt.Rd
-ce8bb515bc6422a640e5a656a9c05a10 *man/uniformmax.Rd
-ee0b70d324af1f4ed599b784ba459a06 *man/unrollPsd.Rd
+5022bab0080b718bae1eab5347e50528 *man/densityMixture.Rd
+ecc6050704152a194433b06485203b09 *man/dtex.kurt.Rd
+6dc937d7338ca1fa399c02e5390528fe *man/fast_ft.Rd
+105117287f46c8b490380a99412d22eb *man/fast_ift.Rd
+3ddf0d0d4e3eb49958d072bd45908f46 *man/fast_mean.Rd
+d1adece5bae596fd71a0a285d7815405 *man/filenameMCMC.Rd
+d860f33ddef48f28ab6859f5055726ce *man/genEpsARC.Rd
+186b02d97acdb0655d5b909d3bcd6ea7 *man/genEpsARMAC.Rd
+27c00a5e205bf9d066d22b9fa3e31d3e *man/genEpsMAC.Rd
+e796ce2c302435a3e37b23f2c6c2d290 *man/generalizedGaussian.alpha.Rd
+ea4a2692e80282e4d71ad60aff5891cc *man/generalizedGaussian.kurtosis.Rd
+2dd6e06472b2a5ce1b0cf8852af792b9 *man/gibbs_AR.Rd
+5a08c739c0343d88e84179d47dd1cbdd *man/gibbs_NP.Rd
+12b9d4d1cb26381da0e0726c5d1fbf1d *man/gibbs_NPC.Rd
+5d294eb05f03f38f1d18f468917faf52 *man/gibbs_pacf.Rd
+d9b9cff216c484f9a8f2015ec853ab26 *man/gibbs_toggle_ar.Rd
+df959d038a31db691167967d14ce750d *man/l_generalizedGaussian.Rd
+2cd014cf13e8d81d2fe9334257090c9f *man/llike.Rd
+0370892c93a380db5fa1e8ee52bd6c03 *man/llike_full.Rd
+6dd944a5bcf08efc6e7713fc721bdc46 *man/llike_pacf.Rd
+3f55eaa780e8ecc6a99ab549f4a9a19d *man/logfuller.Rd
+6e46a23c54462e24e39274dd3360c991 *man/logrosenthal.Rd
+43c1679491f405842c74c8cf4e060dfc *man/lpost.Rd
+ec672b1fb40514c17d053348ab062515 *man/lpost_pacf.Rd
+5263d9bcd4cecfcb9031b5208f13fe28 *man/lprior.Rd
+de40aba0e228feb34203f063d403118c *man/lprior_pacf.Rd
+0f2328c2ea7330a6ded78b6732979b96 *man/mixtureWeight.Rd
+6c6dec110910a5d574ab4ce960c512c9 *man/nll_generalizedGaussian.Rd
+645256a001b8334e28c23a8dc01f36d4 *man/nll_norm.Rd
+5ca2b22a2990d2a17d08569b253e0eb1 *man/nll_norm_unnormalized.Rd
+8b022c13594d267160e6d53d9afa2700 *man/nll_t.Rd
+95a14658fc03dea6ce9aaf2b627b20c5 *man/omegaFreq.Rd
+cc6f4a20438286336dc1e64dff8b1204 *man/pFromV.Rd
+10024533d9e65852fd6a73b65a3177eb *man/pacf2AR.Rd
+fd368aa6764c92c381595bae3d8ff7a6 *man/pacf2ARacv.Rd
+6ab7d4a6d959e4f10566c3087c16138c *man/pacfToAR.Rd
+9ec79ecf0b1c8eddf63e75d129a2b81b *man/plotMCMC.Rd
+8e57521520bacfe7ba6c7ebfe8db6fdf *man/plotPsdEstimate.Rd
+10f3b2135065d6b3dbb4968c3492092e *man/psd_arma.Rd
+4d5c98cb8710479f2c7486d0d4be68bc *man/reduceMemoryStorageMCMC.Rd
+f428995711eefd00ae7e91bcf77b9d40 *man/se_kurt.Rd
+0d98832a5e38fe5a8298e63e9d902c24 *man/uniformmax.Rd
+da4fe5fc3e98576749d9dcc9eb476637 *man/unrollPsd.Rd
 eac368c5d29e3b3ad52f636af556f0ed *src/RcppExports.cpp
-3be4f5584485a98c9d3431ac080fc15b *src/gibbs_core.cpp
+482c201cab2c80b1e577457c9ea33574 *src/gibbs_core.cpp
 f0cc656c96a8b6ef854dddfa660e5c30 *src/gibbs_util.cpp
+7c22b1b500cab25872924e218f1645f5 *src/registerDynamicSymbol.c

NAMESPACE

@@ -31,4 +31,4 @@ importFrom(stats,runif)
 importFrom(stats,sd)
 importFrom(stats,var)
 importFrom(utils,head)
-useDynLib(beyondWhittle)
+useDynLib(beyondWhittle, .registration = TRUE)

R/gibbs_AR.R

@@ -32,7 +32,7 @@
 #' # Use this variable to set the autoregressive model order
 #' ar.order <- 2
 #' 
-#' data <- sqrt(as.numeric(sunspot.year[-length(sunspot.year)]))
+#' data <- sqrt(as.numeric(sunspot.year))
 #' data <- data - mean(data)
 #' 
 #' # If you run the example be aware that this may take several minutes
@@ -91,7 +91,7 @@
 #' @importFrom graphics abline lines plot title
 #' @importFrom stats ar dbeta dgamma dnorm dt fft mad median plot.ts quantile rbeta rgamma rnorm rt runif sd var
 #' @importFrom utils head
-#' @useDynLib beyondWhittle
+#' @useDynLib beyondWhittle, .registration = TRUE
 #' @export
 gibbs_AR <- function(data,
                      Ntotal,

R/gibbs_NP.R

@@ -30,7 +30,7 @@
 #' ## Example 1: Fit the NP model to sunspot data:
 #' ##
 #' 
-#' data <- sqrt(as.numeric(sunspot.year[-length(sunspot.year)]))
+#' data <- sqrt(as.numeric(sunspot.year))
 #' data <- data - mean(data)
 #' 
 #' # If you run the example be aware that this may take several minutes
@@ -83,7 +83,7 @@
 #' cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")
 #' }
 #' @importFrom Rcpp evalCpp
-#' @useDynLib beyondWhittle
+#' @useDynLib beyondWhittle, .registration = TRUE
 #' @export
 gibbs_NP <- function(data,
                      Ntotal,

R/gibbs_NPC.R

@@ -32,13 +32,13 @@
 #' \dontrun{
 #' 
 #' ##
-#' ## Example 1: Fit an AR(p) model to sunspot data:
+#' ## Example 1: Fit a nonparametrically corrected AR(p) model to sunspot data:
 #' ##
 #' 
 #' # Use this variable to set the autoregressive model order
 #' ar.order <- 2
 #' 
-#' data <- sqrt(as.numeric(sunspot.year[-length(sunspot.year)]))
+#' data <- sqrt(as.numeric(sunspot.year))
 #' data <- data - mean(data)
 #' 
 #' # If you run the example be aware that this may take several minutes
@@ -61,7 +61,7 @@
 #' 
 #' 
 #' ##
-#' ## Example 2: Fit an AR(p) model to high-peaked AR(1) data
+#' ## Example 2: Fit a nonparametrically corrected AR(p) model to high-peaked AR(1) data
 #' ##
 #' 
 #' # Use this variable to set the autoregressive model order
@@ -96,7 +96,7 @@
 #' cat("IAE=", mean(abs(mcmc$psd.median-psd_true)[-1]) , sep="")
 #' }
 #' @importFrom Rcpp evalCpp
-#' @useDynLib beyondWhittle
+#' @useDynLib beyondWhittle, .registration = TRUE
 #' @export
 gibbs_NPC <- function(data,
                       Ntotal,

R/intern__gibbs_core.R

@@ -88,59 +88,66 @@ arma_conditional <- function(zt, ar, ma, nll_fun, sigma2, ...) {
 #' Log corrected parametric likelihood
 #' @keywords internal
 llike <- function(omega, FZ, ar, ma, v, w, k, tau, corrected, toggle.q, pdgrm, db.list, dist, nll_fun, ex.kurt, beta, f.alpha) {
-
+  
   # Calculates Whittle or corrected log-likelihood (assume n even) for Gaussian errors
   n <- length(FZ)
-
+  
   # Un-normalised PSD (defined on [0, 1])
   q.psd <- qpsd(omega, v, w, k, db.list)$psd
   q <- unrollPsd(q.psd, n)
-
+  
   # Normalised PSD (defined on [0, pi])
   f <- tau * q
   #f_param <- psd_arma(pi*omega, ar, ma,sigma2.ex)
-
-
+
+  # do not consider the following frequencies in the likelihood,
+  # as they correspond to the mean (or alternating mean)
+  if (n %% 2) {
+    excludedFrequecies <- c(1,n)
+  } else {
+    excludedFrequecies <- 1
+  }
+
   # Corrected log-likelihood
   if (corrected == TRUE) {
-
-#     if (dist=="student" && ex.kurt != 0) {
-#       df <- 6 / ex.kurt + 4
-#       sigma2 <- df / (df - 2) # sigma2 of likelihood (!=1 only possible for t data, have to correct it here)
-#     } else {
-#       sigma2 <- 1
-#     }
+    
+    #     if (dist=="student" && ex.kurt != 0) {
+    #       df <- 6 / ex.kurt + 4
+    #       sigma2 <- df / (df - 2) # sigma2 of likelihood (!=1 only possible for t data, have to correct it here)
+    #     } else {
+    #       sigma2 <- 1
+    #     }
     sigma2 <- 1
-
+    
     if (toggle.q) {
-      C <- sqrt(unrollPsd(psd_arma(pi*omega,ar,ma,1),n)^(1-f.alpha) / f)
+      C <- sqrt(unrollPsd(psd_arma(pi*omega,ar,ma,1),n)^(1-f.alpha) / f) 
     } else {
       C <- sqrt(unrollPsd(psd_arma(pi*omega,ar,ma,1),n) / f) #Cn(n, f)
     }
-
+    
     # Input for ARMA parametric likelihood - Inverse FT
     FCFZ <- fast_ift(C * FZ)
-
+    
     # Calculate ARMA parametric conditional likelihood
     if (dist == "generalized") {
       ex.kurt <- beta # QUICK HACK: Abuse ex.kurt to parse beta to nll_fun
     }
     p.arma <- arma_conditional(FCFZ, ar, ma, nll_fun, sigma2, ex.kurt)
-
-    # Corrected Gaussian log-likelihood
+    
+    # Corrected Gaussian log-likelihood 
     # llike <- sum(log(C[-c(1,n)])) - p.arma   # Note: The minus sign here.
-    llike <- sum(log(C[-c(1,n)])) - p.arma   # Note: The minus sign here.
-
+    llike <- sum(log(C[-excludedFrequecies])) - p.arma   # Note: The minus sign here.
+    
   }
-
+  
   # Whittle log-likelihood
   if (corrected == FALSE) {
-    llike <- -sum(log(f[2:(n - 1)] * 2 * pi) + pdgrm[2:(n - 1)] / (f[2:(n - 1)] * 2 * pi))  # Remove first and last here
+    llike <- -sum(log(f[-excludedFrequecies] * 2 * pi) + pdgrm[-excludedFrequecies] / (f[-excludedFrequecies] * 2 * pi))
     llike <- llike / 2
-  }
-
+  } 
+  
   return(llike)
-
+  
 }
 
 #' Log corrected posterior
@@ -151,11 +158,55 @@ lpost <- function(omega, FZ, ar, ma, v, w, k, tau,
                   dist, nll_fun, ex.kurt, kurt.lambda, beta, beta.alpha, beta.beta, f.alpha, rho, rho.alpha, rho.beta) {
 
   # Unnormalised log posterior
-  lp <- llike(omega, FZ, ar, ma, v, w, k, tau, corrected, toggle.q, pdgrm, db.list, dist, nll_fun, ex.kurt, beta, f.alpha) +
-    lprior(v, w, k, tau, M, g0.alpha, g0.beta, k.theta, tau.alpha, tau.beta, dist, ex.kurt, kurt.lambda, beta, beta.alpha, beta.beta, f.alpha, corrected, rho, rho.alpha, rho.beta)
-
-  return(lp)
-
+  ll <- llike(omega, FZ, ar, ma, v, w, k, tau, corrected, toggle.q, pdgrm, db.list, dist, nll_fun, ex.kurt, beta, f.alpha)
+  lp <- lprior(v, w, k, tau, M, g0.alpha, g0.beta, k.theta, tau.alpha, tau.beta, dist, ex.kurt, kurt.lambda, beta, beta.alpha, beta.beta, f.alpha, corrected, rho, rho.alpha, rho.beta)
+
+  if (is.na(ll)) {
+    ll_params <- list(omega=omega, 
+                      FZ=FZ, 
+                      ar=ar, 
+                      ma=ma, 
+                      v=v, 
+                      w=w, 
+                      k=k,
+                      tau=tau, 
+                      corrected=corrected, 
+                      toggle.q=toggle.q, 
+                      pdgrm=pdgrm, 
+                      dist=dist, 
+                      nll_fun=nll_fun, 
+                      ex.kurt=ex.kurt, 
+                      beta=beta, 
+                      f.alpha=f.alpha)
+    TXT <- paste("Likelihood evaluated as NA. Parameters: ", ll_params)
+    stop(TXT)
+  }
+  if (is.na(lp)) {
+    lp_params <- list(v=v, 
+                      w=w, 
+                      k=k, 
+                      tau=tau, 
+                      M=M, 
+                      g0.alpha=g0.alpha, 
+                      g0.beta=g0.beta, 
+                      k.theta=k.theta, 
+                      tau.alpha=tau.alpha, 
+                      tau.beta=tau.beta, 
+                      dist=dist, 
+                      ex.kurt=ex.kurt, 
+                      kurt.lambda=kurt.lambda, 
+                      beta=beta, 
+                      beta.alpha=beta.alpha, 
+                      beta.beta=beta.beta, 
+                      f.alpha=f.alpha, 
+                      corrected=corrected, 
+                      rho=rho, 
+                      rho.alpha=rho.alpha, 
+                      rho.beta=rho.beta)
+    TXT <- paste("Prior evaluated as NA. Parameters: ", lp_params)
+    stop(TXT)
+  }
+  return(ll+lp)
 }
 
 #' Log prior of Bernstein-Dirichlet mixture and parametric working model -- all unnormalized