#' This function takes into account a posterior pilot sample to calibrate
#' the proposals for the weights
#'
#' @importFrom Rcpp evalCpp
#' @importFrom expm sqrtm
#' @importFrom stats median rnorm runif rgamma
#' @importFrom fda create.bspline.basis bsplinepen
#' @useDynLib psplinePsd, .registration = TRUE
#' @keywords internal
gibbs_pspline_postProposal <- function(data,
Ntotal,
burnin,
thin = 1,
tau.alpha, # specified in psd
tau.beta, # specified in psd
phi.alpha, # specified in psd
phi.beta, # specified in psd
delta.alpha,# specified in psd
delta.beta, # specified in psd
k, # specified in psd
eqSpacedKnots, # specified in psd
degree, # specified in psd
diffMatrixOrder, # specified in psd
printIter = 100,
psd,
add = FALSE) {
### extracting psd information ###
if(is.null(psd)){
stop("Include an output from gibbs_psline function");
}else{
if(class(psd) != "psd")stop("'psd' must be a psd object");
cat("Posterior samples used to calibrate the proposals for the weights", "\n")
}
k = psd$anSpecif$k;
degree = psd$anSpecif$degree;
eqSpacedKnots = psd$anSpecif$eqSpacedKnots;
diffMatrixOrder = psd$anSpecif$diffMatrixOrder;
tau.alpha = psd$anSpecif$tau.alpha;
tau.beta = psd$anSpecif$tau.beta;
phi.alpha = psd$anSpecif$phi.alpha;
phi.beta = psd$anSpecif$phi.beta;
delta.alpha = psd$anSpecif$delta.alpha;
delta.beta = psd$anSpecif$delta.beta;
if (burnin >= Ntotal) stop("burnin must be less than Ntotal")
if (any(c(Ntotal, thin) %% 1 != 0) || any(c(Ntotal, thin) <= 0)) stop("Ntotal and thin must be strictly positive integers")
if ((burnin %% 1 != 0) || (burnin < 0)) stop("burnin must be a non-negative integer")
if (any(c(tau.alpha, tau.beta) <= 0)) stop("tau.alpha and tau.beta must be strictly positive")
if (any(c(phi.alpha, phi.beta) <= 0)) stop("phi.alpha and phi.beta must be strictly positive")
if (any(c(delta.alpha, delta.beta) <= 0)) stop("delta.alpha and delta.beta must be strictly positive")
if(!is.logical(add))stop("add must be a logical value");
cat(paste("Number of B-splines k=", k, sep=""), "\n");
if( (Ntotal - burnin)/thin < k){
stop("Change specifications: Either increase Ntotal or decrease thin parameters, otherwise the covariance matrix of the weights is singular");
}
### ###
n <- length(data);
# Which boundary frequencies to remove from likelihood computation and tau sample
if (n %% 2) { # Odd length time series
bFreq <- 1 # Remove first
}
else { # Even length time series
bFreq <- c(1, n) # Remove first and last
}
if( (printIter<=0) || (printIter %% 1 != 0) )stop("printIter must be a positive integer value");
if(!(k-2 >= diffMatrixOrder))stop("diffMatrixOrder must be lower than or equal to k-2");
# Tolerance for mean centering
tol <- 1e-4;
# Mean center
if (abs(mean(data)) > tol) {
data <- data - mean(data)
warning("Data has been mean-centred")
}
# Optimal rescaling to prevent numerical issues
rescale = stats::sd(data);
data = data / rescale; # Data now has standard deviation 1
FZ <- fast_ft(data); # FFT data to frequency domain. NOTE: Must be mean-centred.
pdgrm <- abs(FZ) ^ 2; # Periodogram: NOTE: the length is n here.
omega <- 2 * (1:(n / 2 + 1) - 1) / n; # Frequencies on unit interval
lambda <- pi * omega; # Angular frequencies on [0, pi]
# number of samples to be stored
N = round(Ntotal/thin);
# Open objects for storage
tau <- rep(NA, N);
phi <- rep(NA, N); # scalar in Multivariate Normal
delta <- rep(NA, N);
######################################
### Recycling gibbs_pspline output ###
######################################
# Proposal calibration
muV = apply(psd$V, 1, mean); # mean vector - posterior samples
covV = stats::cov(t(psd$V)); # covariance matrix - posterior samples
sqrt.covV = expm::sqrtm(covV);
# Starting values
if( add ==TRUE){
pilot_N = length(psd$tau);
tau[1] <- psd$tau[pilot_N];
delta[1] <- psd$delta[pilot_N];
phi[1] <- psd$phi[pilot_N];
v <- psd$V[, pilot_N];
}else{
tau[1] <- stats::median(psd$tau); # stats::var(data) / (2 * pi);
delta[1] <- stats::median(psd$delta);# delta.alpha / delta.beta;
phi[1] <- stats::median(psd$phi); # phi.alpha/(phi.beta * delta[1]);
v <- apply(psd$V, 1, median);
}
v <- solve(sqrt.covV) %*% (v - muV); # detransforming
V <- matrix(NA, nrow = length(v), ncol = N);
V[, 1] <- v;
#V = matrix(v, ncol = 1); # To store beta values, see page 6 paper
# Fixed knot number & location => a single calculation of B-spline densities
knots <- knotLoc(data = data, k = k, degree = degree, eqSpaced = eqSpacedKnots);
db.list <- dbspline(omega, knots, degree);
# Difference Matrix
if(eqSpacedKnots == TRUE){
P <- diffMatrix(k-1, d = diffMatrixOrder); # Third order penalty
P <- t(P) %*% P;
}else{
if(degree <= diffMatrixOrder){
stop("bsplinepen function: penalty order must be lower than the bspline density degree");
}
nknots <- length(knots);
basisobj <- fda::create.bspline.basis(c(0, knots[nknots-1]), norder = degree + 1,
nbasis = k - 1, breaks = knots[-nknots]);
P <- fda::bsplinepen(basisobj, Lfdobj = diffMatrixOrder);
P <- P/norm(P);
}
epsilon <- 1e-6; #1e-6;
P <- P + epsilon * diag(dim(P)[2]); # P^(-1)=Sigma (Covariance matrix)
# Store log likelihood
ll.trace <- NULL;
Count <- NULL; # acceptance probability
sigma <- 1; # variance of proposal distb for weights
count <- 0.4; # starting value for acc pbb - optimal value
k1 <- k - 1;
# Random values
Zs <- stats::rnorm((Ntotal-1)*k1);
Zs <- matrix(Zs, nrow = Ntotal-1, ncol = k1);
Us <- log(stats::runif((Ntotal-1)*k1, 0, 1));
Us <- matrix(Us, nrow = Ntotal-1, ncol = k1);
# Initial values for the proposals
phi.store <- phi[1];
tau.store <- tau[1];
delta.store <- delta[1];
V.store <- V[, 1];
ptime <- proc.time()[1]
# Metropolis-within-Gibbs sampler
for(j in 1:(N-1)){
adj <- (j - 1) * thin;
V.star <- V.store; # proposal value
aux <- sample(k1); # positions to be changed in the thining loop
# Thinning
for (i in 1:thin) {
iter <- i + adj;
if (iter %% printIter == 0) {
cat(paste("Iteration", iter, ",", "Time elapsed",
round(as.numeric(proc.time()[1] - ptime) / 60, 2),
"minutes"), "\n")
}
v.store <- sqrt.covV %*% V.store + muV;
f.store <- lpost(omega,
FZ,
k,
v.store, # parameter
tau.store, # parameter
tau.alpha,
tau.beta,
phi.store, # parameter
phi.alpha,
phi.beta,
delta.store, # parameter
delta.alpha,
delta.beta,
P,
pdgrm,
degree,
db.list)
##############
### WEIGHT ###
##############
#V.star = V.store;
#aux = sample(k1);
# tunning proposal distribution
if(count < 0.30){ # increasing acceptance pbb
sigma <- sigma * 0.90; # decreasing proposal moves
}else if(count > 0.50){ # decreasing acceptance pbb
sigma <- sigma * 1.1; # increasing proposal moves
}
count <- 0; # acceptance probability
for(g in 1:k1){
pos <- aux[g];
#sigma <- Sigma[pos]; # Another approach, individual variance
V.star[pos] <- V.store[pos] + sigma * Zs[iter, g];
v.star = sqrt.covV %*% V.star + muV;
f.V.star <- lpost(omega,
FZ,
k,
v.star, # proposal value
tau.store,
tau.alpha,
tau.beta,
phi.store,
phi.alpha,
phi.beta,
delta.store,
delta.alpha,
delta.beta,
P,
pdgrm,
degree,
db.list);
#Accept/reject
alpha1 <- min(0, f.V.star - f.store); # log acceptance ratio
if(Us[iter, g] < alpha1) {
V.store[pos] <- V.star[pos]; # Accept V.star
f.store <- f.V.star;
count <- count + 1; # ACCEPTANCE PROBABILITY
v.store <- v.star;
}else {
V.star[pos] = V.store[pos]; # reseting proposal value
}
} # End updating weights
count <- count / k1;
Count[iter] <- count; # acceptance probability
###########
### phi ###
###########
#v = sqrt.covV %*% V.store + muV;
v <- v.store;
phi.store <- stats::rgamma(1, shape = k1/2 + phi.alpha,
rate = phi.beta * delta.store + t(v) %*% P %*% v / 2);
#############
### delta ###
#############
delta.store <- stats::rgamma(1, shape = phi.alpha + delta.alpha,
rate = phi.beta * phi.store + delta.beta);
###########
### tau ###
###########
# Directly sample tau from conjugate Inverse-Gamma density
q.psd <- qpsd(omega, k, v.store, degree, db.list);
q <- unrollPsd(q.psd, n)
# Note: (n - 1) and (n - 2) here. Remove the first and last terms for even and first for odd
if (n %% 2) { # Odd length series
tau.store <- 1 / stats::rgamma(1, tau.alpha + (n - 1) / 2,
tau.beta + sum(pdgrm[-bFreq] / q[-bFreq]) / (2 * pi) / 2)
}
else{ # Even length series
tau.store <- 1 / stats::rgamma(1, tau.alpha + (n - 2) / 2,
tau.beta + sum(pdgrm[-bFreq] / q[-bFreq]) / (2 * pi) / 2)
}
} # END: thining
######################
### Storing values ###
######################
phi[j+1] <- phi.store;
delta[j+1] <- delta.store;
tau[j+1] <- tau.store;
V[, j+1] <- V.store;
} # END: MCMC loop
# Discarding burn-in period
keep <- seq(round(burnin/thin) + 1, N);
tau <- tau[keep];
phi <- phi[keep];
delta<- delta[keep];
V <- V[, keep];
# converting to actual V
V = apply(V,2, function(x) sqrt.covV %*% x + muV);
##############################
### Compute log likelihood ###
##############################
for(i in 1:length(keep)){
ll.trace <- c(ll.trace,
llike(omega, FZ, k, V[,i], tau[i], pdgrm, degree, db.list));
}
fpsd.sample <- log.fpsd.sample <- matrix(NA, nrow = length(omega), ncol = length(keep));
# Store PSDs
for (isample in 1:length(keep)) {
q.psd <- qpsd(omega, k, V[, isample], degree, db.list);
fpsd.sample[, isample] <- tau[isample] * q.psd;
log.fpsd.sample[, isample] <- logfuller(fpsd.sample[, isample]); # Create transformed version
}
if(add == TRUE){
cat("Pilot and current analyses have been merged", "\n");
count <- c(psd$count * k, Count); # k factor is cancelled return output
# db.list is equal in both analyses
delta <- c(psd$delta, delta);
aux <- psd$fpsd.sample / (rescale^2);
fpsd.sample <- cbind(aux, fpsd.sample);
log.fpsd.sample <- cbind(apply(aux, 2, logfuller), log.fpsd.sample);
ll.trace <- c(psd$ll.trace, ll.trace);
# pdgrm is equal in both analyses
phi <- c(psd$phi, phi);
# psd.mean, psd.median, psd.p05, psd.p95, psd.u95
# are calculated below from fpsd.sample
tau <- c(psd$tau, tau);
V <- cbind(psd$V, V);
}
# Compute point estimates and 90% Pointwise CIs
psd.median <- apply(fpsd.sample, 1, stats::median);
psd.mean <- apply(fpsd.sample, 1, mean);
psd.p05 <- apply(fpsd.sample, 1, stats::quantile, probs=0.05);
psd.p95 <- apply(fpsd.sample, 1, stats::quantile, probs=0.95);
# Transformed versions of these for uniform CI construction
log.fpsd.s <- apply(log.fpsd.sample, 1, stats::median);
log.fpsd.mad <- apply(log.fpsd.sample, 1, stats::mad);
log.fpsd.help <- apply(log.fpsd.sample, 1, uniformmax);
log.Cvalue <- stats::quantile(log.fpsd.help, 0.9);
# Compute Uniform CIs
psd.u95 <- exp(log.fpsd.s + log.Cvalue * log.fpsd.mad);
psd.u05 <- exp(log.fpsd.s - log.Cvalue * log.fpsd.mad);
###########
### DIC ###
###########
tau_mean <- mean(tau);
v_means <- unname(apply(V, 1, mean)); # V was converted above
l <- llike(omega, FZ, k, v = v_means,
tau = tau_mean, pdgrm, degree, db.list);
ls <- apply(rbind(tau, V), 2, function(x){
llike(omega, FZ, k, v = x[-1],
tau = x[1], pdgrm, degree, db.list)});
ls <- unname(ls);
# http://kylehardman.com/BlogPosts/View/6
# DIC = -2 * (l - (2 * (l - mean(ls))));
D_PostMean <- -2 * l;
D_bar <- -2 * mean(ls);
pd <- D_bar - D_PostMean;
DIC <- list(DIC = 2 * D_bar - D_PostMean, pd = pd);
# Compute periodogram
N <- length(psd.median); # N = (n + 1) / 2 (ODD) or N = n / 2 + 1 (EVEN)
pdgrm <- (abs(stats::fft(data)) ^ 2 / (2 * pi * n))[1:N];
# storing analysis specifications
anSpecif <- list(k = k, n = n, degree = degree, FZ = FZ,
eqSpacedKnots = eqSpacedKnots, diffMatrixOrder = diffMatrixOrder,
tau.alpha = tau.alpha, tau.beta = tau.beta,
phi.alpha = phi.alpha, phi.beta = phi.beta,
delta.alpha = delta.alpha, delta.beta = delta.beta);
# List to output
output <- list(psd.median = psd.median * rescale ^ 2,
psd.mean = psd.mean * rescale ^ 2,
psd.p05 = psd.p05 * rescale ^ 2,
psd.p95 = psd.p95 * rescale ^ 2,
psd.u05 = psd.u05 * rescale ^ 2,
psd.u95 = psd.u95 * rescale ^ 2,
fpsd.sample = fpsd.sample * rescale ^ 2,
anSpecif = anSpecif,
n = n,
tau = tau,
phi = phi,
delta = delta,
V = V,
ll.trace = ll.trace,
pdgrm = pdgrm * rescale ^ 2,
db.list = db.list,
DIC = DIC,
count = Count); # Acceptance probability
class(output) <- "psd" # Assign S3 class to object
return(output); # Return output
} # Close function
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.