R/mcmc_update.R

In PPgp/bayesTFR: Bayesian Fertility Projection

mcmc.update.abS <- function(what, eps_Tc_temp, mcmc) {
  # 'what' is one of ('a', 'b', 'S')
  var.index <- (1:3)[what == c('a', 'b', 'S')]
  abS.values <- list(a=mcmc$a_sd, b=mcmc$b_sd, S=mcmc$S_sd)
  var.value <- abS.values[[var.index]]
  var.low <- mcmc$meta[[paste(what, 'low', sep='.')]]
  var.up <- mcmc$meta[[paste(what, 'up', sep='.')]]
  var.width <- mcmc$meta[[paste(what, 'width', sep='.')]]
  var.name <- paste(what, 'sd', sep='_')
  
  ##############################################################
  # slice sampling starts with getting the z = log_cond full likelihood - exp(1)
  # and then finding a new estimate such that its log_cond_full_likelihood is higher than z
  # the new estimate is found within interval around current estimate
  ##############################################################
  
  z <- log_cond_abf_sd(mcmc$add_to_sd_Tc, mcmc$const_sd, mcmc$sigma0, eps_Tc_temp, 
                       mcmc$const_sd_dummie_Tc, mcmc$meta$sigma0.min) - rexp(1)
  # if (what == 'a') cat(mcmc$finished.iter, ":", z, "a:", mcmc$a_sd, "\n")
  v <- runif(1)   
  interval <- c(max(var.value - v*var.width, var.low),
                min(var.value + (1-v)*var.width,var.up))
  
  add_to_sd_Tc_prop <- matrix(NA, nrow(mcmc$add_to_sd_Tc), ncol(mcmc$add_to_sd_Tc))
  for (country in 1:mcmc$meta$nr_countries){
    add_to_sd_Tc_prop[1:length(mcmc$data.list[[country]]), country] <- (mcmc$data.list[[country]] - abS.values$S)*
      ifelse(mcmc$data.list[[country]] > abS.values$S, 
             -abS.values$a, abS.values$b)
  }  
  #while (TRUE){
  for(i in 1:50) {
    var_prop <- runif(1,interval[1], interval[2])
    abS.values[[what]] <- var_prop
    for (country in 1:mcmc$meta$nr_countries){
      add_to_sd_Tc_prop[1:length(mcmc$data.list[[country]]), country] <- (mcmc$data.list[[country]] - abS.values$S)*
        ifelse(mcmc$data.list[[country]] > abS.values$S, 
               -abS.values$a, abS.values$b)
    }  
    like <- log_cond_abf_sd(add_to_sd_Tc_prop, mcmc$const_sd, mcmc$sigma0, eps_Tc_temp, 
                            mcmc$const_sd_dummie_Tc, mcmc$meta$sigma0.min)     
    if (like >= z) {
      mcmc[[var.name]] <- var_prop
      mcmc$add_to_sd_Tc <- add_to_sd_Tc_prop
      return()
    } else {
      # shrink interval
      if (var_prop < var.value){
        interval[1] <- var_prop
      } else {
        interval[2] <- var_prop
      }
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else
  }
  mcmc[[var.name]] <- var_prop
  mcmc$add_to_sd_Tc <- add_to_sd_Tc_prop
  warning("Cannot find likelihood increase for ", what, ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate. = TRUE)
}

mcmc.update.sigma0const <- function(what, log.like.func, eps_Tc_temp, mcmc) {
  # 'what' is one of ('sigma0', 'const')
  var.index <- (1:2)[what == c('sigma0', 'const')]
  var.values <- list(sigma0=mcmc$sigma0, const=mcmc$const_sd)
  var.name <- c('sigma0', 'const_sd')[var.index]
  var.value <- var.values[[var.index]]
  var.low <- mcmc$meta[[paste(what, 'low', sep='.')]]
  var.up <- mcmc$meta[[paste(what, 'up', sep='.')]]
  var.width <- mcmc$meta[[paste(what, 'width', sep='.')]]
  
  z <- eval(call(log.like.func, mcmc$add_to_sd_Tc, var.values[['const']], var.values[['sigma0']], 
                 eps_Tc_temp, mcmc$const_sd_dummie_Tc, mcmc$meta$sigma0.min)) - rexp(1)
  v <- runif(1)
  interval <- c(max(var.value - v*var.width, var.low),
                min(var.value + (1-v)*var.width,var.up))
  
  #while (TRUE){
  for(i in 1:50) {
    var_prop <- runif(1,interval[1], interval[2])
    var.values[[what]] <- var_prop
    like <- eval(call(log.like.func, mcmc$add_to_sd_Tc, var.values[['const']], var.values[['sigma0']], 
                      eps_Tc_temp, mcmc$const_sd_dummie_Tc, mcmc$meta$sigma0.min))
    if (like >= z) {
      mcmc[[var.name]] <- var_prop
      return()
    } else {
      # shrink interval
      if (var_prop < var.value){
        interval[1] <- var_prop
      } else {
        interval[2] <- var_prop
      }
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else
  }
  mcmc[[var.name]] <- var_prop
  warning("Cannot find likelihood increase for ", what, ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate.=TRUE)
}

mcmc.update.rho.phase2 <- function(mcmc, matrix.name) {
  var.value <- mcmc$rho.phase2
  var.low <- 0
  var.up <- 0.9
  var.width <- 0.1
  
  z <- sum(dnorm(mcmc$eps_Tc, mean = 0, sd =  mcmc$sd_Tc, log = TRUE), na.rm = TRUE) - rexp(1)
  v <- runif(1)
  interval <- c(max(var.value - v*var.width, var.low), min(var.value + (1-v)*var.width,var.up))
  
  #while (TRUE){
  for(i in 1:50) {
    var_prop <- runif(1,interval[1], interval[2])
    eps_prop <- get_eps_T_all(mcmc, matrix.name=matrix.name, rho.phase2=var_prop)
    like <- sum(dnorm(eps_prop, mean = 0, sd =  mcmc$sd_Tc, log = TRUE), na.rm = TRUE)
    if (like >= z) {
      mcmc$rho.phase2 <- var_prop
      mcmc$eps_Tc <- eps_prop
      return()
    } else {
      # shrink interval
      if (var_prop < var.value){
        interval[1] <- var_prop
      } else {
        interval[2] <- var_prop
      }
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else
  }
  mcmc$rho.phase2 <- var_prop
  mcmc$eps_Tc <- eps_prop
  warning("Cannot find likelihood increase for ", "rho phase2", ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate. = TRUE)
}

mcmc.update.abSsigma0const <- function(mcmc, id.not.early.index) {
  # updates a_sd, b_sd, f_sd, sigma0, and const_sd
  #################################
  
  # propose a new a, b, f sigma0 and c (each separately)
  # then calculate the new distortions to get its log_cond_like
  # then check if you want to accept the new parameter, or shrink the interval and repeat
  # if accepted, update the parameters itself and add_to_sd_Tc matrix
  # and at the end, update sd_Tc (as this is needed when updating other parameters)
  
  # temp put NAs in eps_ctau's
  eps_Tc_temp = mcmc$eps_Tc
  # put NAs for countries that are not included in the estimation
  eps_Tc_temp[(mcmc$meta$nr_countries_estimation+1):mcmc$meta$nr_countries] <- NA
  
  # put NAs on tau_c spots for id_not_early countries
  eps_Tc_temp[id.not.early.index] <- NA
  
  mcmc.update.abS('a', eps_Tc_temp, mcmc)
  mcmc.update.abS('b', eps_Tc_temp, mcmc)
  mcmc.update.abS('S', eps_Tc_temp, mcmc)
  mcmc.update.sigma0const('sigma0', 'log_cond_sigma0', eps_Tc_temp, mcmc)
  mcmc.update.sigma0const('const', 'log_cond_const_sd', eps_Tc_temp, mcmc)
  mcmc$sd_Tc <- ifelse(mcmc$const_sd_dummie_Tc==1, mcmc$const_sd, 1)*
    ifelse((mcmc$sigma0 + mcmc$add_to_sd_Tc)>0, mcmc$sigma0 + mcmc$add_to_sd_Tc, 
           mcmc$meta$sigma0.min)
}

##############################################################################
mcmc.update.Triangle_c4 <- function(country, mcmc, ...) {
  # within country loop, updates Triangle_c4
  # thus also the eps
  # (thus also lambda_c and the NAs in the eps!)
  Triangle_c4_trans <- log(max(mcmc$Triangle_c4[country] - mcmc$meta$Triangle_c4.low, 1e-20)/
                             max(mcmc$meta$Triangle_c4.up - mcmc$Triangle_c4[country], 1e-20))

  epsT.idx <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country]-1)
  # exclude indices with extreme epsT.idx
  raw.outliers <- mcmc$meta$indices.outliers[[as.character(country)]]
  if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2) 
    raw.outliers <- sort(unique(c(raw.outliers, raw.outliers+1)))
  epsT.idx <- setdiff(epsT.idx, raw.outliers)
  
  lepsT.idx <- length(epsT.idx)
  #        z <- (log_cond_Triangle_c4_trans(Triangle_c4_trans, 
  #                mcmc$eps_Tc[epsT.idx,country],
  #                mcmc$sd_Tc[epsT.idx,country],
  #                mcmc$mean_eps_Tc[epsT.idx,country], mcmc$Triangle4, mcmc$delta4) - rexp(1))
  log_cond <- 0.0
  logcondt <- .C("log_cond_Triangle_c4_trans", Triangle_c4_trans, 
                 mcmc$eps_Tc[epsT.idx,country],
                 mcmc$sd_Tc[epsT.idx,country],
                 mcmc$mean_eps_Tc[epsT.idx,country],
                 lepsT.idx, mcmc$Triangle4, mcmc$delta4, log_cond=log_cond,
                 PACKAGE = "bayesTFR")
  z <-  logcondt$log_cond - rexp(1)
  #		stop('')
  v <- runif(1)
  
  interval <- c(Triangle_c4_trans - v*mcmc$meta$Triangle_c4.trans.width, 
                Triangle_c4_trans + (1-v)*mcmc$meta$Triangle_c4.trans.width)
  
  #while (TRUE){
  for(i in 1:50) {
    Triangle_c4_trans_prop <- runif(1,interval[1], interval[2])
    Triangle_c4_prop <- (mcmc$meta$Triangle_c4.up*exp(Triangle_c4_trans_prop) +mcmc$meta$Triangle_c4.low)/
      (1+exp(Triangle_c4_trans_prop))
    theta_prop <-  c((mcmc$U_c[country]-Triangle_c4_prop)*
                       exp(mcmc$gamma_ci[country,])/
                       sum(exp(mcmc$gamma_ci[country,])), 
                     Triangle_c4_prop, 
                     mcmc$d_c[country])
    eps_T_prop <- get.eps.T(theta_prop,country, mcmc$meta, ...)
    like <- .C("log_cond_Triangle_c4_trans", Triangle_c4_trans_prop, eps_T_prop,
               mcmc$sd_Tc[epsT.idx,country],
               mcmc$mean_eps_Tc[epsT.idx,country], lepsT.idx, mcmc$Triangle4, mcmc$delta4,
               log_cond=log_cond, PACKAGE = "bayesTFR")$log_cond
    if (like >= z) {
      mcmc$eps_Tc[epsT.idx, country] <- eps_T_prop
      mcmc$Triangle_c4[country] <- Triangle_c4_prop
      return()
    } else {
      # shrink interval
      if (Triangle_c4_prop < mcmc$Triangle_c4[country]){
        interval[1] <- Triangle_c4_trans_prop
      } else {
        interval[2] <- Triangle_c4_trans_prop
      }
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else
  }
  mcmc$eps_Tc[epsT.idx, country] <- eps_T_prop
  mcmc$Triangle_c4[country] <- Triangle_c4_prop
  warning("Cannot find likelihood increase for Triangle_c4_", mcmc$meta$regions$country_code[country], ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate. = TRUE)
}



############################################
mcmc.update.gamma <- function(country, mcmc, ...) {
  # within country loop, update gamma_c's
  #################################
  # MH-step: propose new vector of gamma's, 
  # calculate the new distortions and log-full-cond-likelihood
  # accept or reject, if accepted update gamma's and the distortions
  
  #################################
  # block update, propose new gamma_ci
  #################################
  gamma_prop <- rmvnorm(1, mcmc$gamma_ci[country,], mcmc$meta$proposal_cov_gammas_cii[country,,])
  pci_prob <- exp(gamma_prop)/sum(exp(gamma_prop))
  theta_prop <- c(pci_prob*(mcmc$U_c[country] - mcmc$Triangle_c4[country]), 
                  mcmc$Triangle_c4[country], mcmc$d_c[country]) 
  eps_T_prop <- get.eps.T(theta_prop, country, mcmc$meta, ...)

  idx <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country]-1)
  # exclude indices with extreme eps in idx
  raw.outliers <- mcmc$meta$indices.outliers[[as.character(country)]]
  if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2) 
    raw.outliers <- sort(unique(c(raw.outliers, raw.outliers+1)))
  
  idx <- setdiff(idx, raw.outliers)
  
  prob_accept <- exp(log_like_gammas(gamma_prop,eps_T_prop,
                                     mcmc$sd_Tc[idx, country], mcmc$mean_eps_Tc[idx, country], 
                                     mcmc$alpha, mcmc$delta) - 
                       log_like_gammas(mcmc$gamma_ci[country,,drop=FALSE],
                                       mcmc$eps_Tc[idx, country], mcmc$sd_Tc[idx, country], 
                                       mcmc$mean_eps_Tc[idx, country], 
                                       mcmc$alpha, mcmc$delta) )
  
  if (runif(1) < prob_accept){
    mcmc$gamma_ci[country,] <- gamma_prop
    mcmc$eps_Tc[idx, country] <- eps_T_prop
  }
  return()
}



mcmc.update.d <- function(country, mcmc, ...) {
  # if accepted, update d_c and the distortions
  d_trans <- log((mcmc$d_c[country] - mcmc$meta$d.low)/(mcmc$meta$d.up - mcmc$d_c[country]))
  idx <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country]-1)
  # Exclude extreme indices in computing loglik
  raw.outliers <- mcmc$meta$indices.outliers[[as.character(country)]]
  if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2) 
    raw.outliers <- sort(unique(c(raw.outliers, raw.outliers+1)))
  
  idx <- setdiff(idx, raw.outliers)
  z <- (log_cond_d_trans(d_trans, 
                         mcmc$eps_Tc[idx, country],
                         mcmc$sd_Tc[idx, country],
                         mcmc$mean_eps_Tc[idx, country],
                         mcmc$psi, mcmc$chi)
        - rexp(1))
  
  v <- runif(1)
  interval <- c(d_trans - v*mcmc$meta$d.trans.width, 
                d_trans + (1-v)*mcmc$meta$d.trans.width)
  
  theta_prop <-  c((mcmc$U_c[country]-mcmc$Triangle_c4[country])*
                     exp(mcmc$gamma_ci[country,])/
                     sum(exp(mcmc$gamma_ci[country,])), mcmc$Triangle_c4[country], 
                   mcmc$d_c[country])
  
  #while (TRUE){
  for(i in 1:50) {
    d_trans_prop <- runif(1,interval[1], interval[2])
    d_prop <- (mcmc$meta$d.up*exp(d_trans_prop) +mcmc$meta$d.low)/(1+exp(d_trans_prop))
    eps_T_prop <- get.eps.T(c(theta_prop[-5], d_prop),country, mcmc$meta, ...)
    if ((like <- log_cond_d_trans(d_trans_prop,eps_T_prop, 
                         mcmc$sd_Tc[idx, country],
                         mcmc$mean_eps_Tc[idx, country],
                         mcmc$psi, mcmc$chi)) >= z) {
      mcmc$eps_Tc[idx, country] <- eps_T_prop
      mcmc$d_c[country] <- d_prop
      return()
    } else {
      # shrink interval
      if (d_prop < mcmc$d_c[country])
        interval[1] <- d_trans_prop
      else 
        interval[2] <- d_trans_prop
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else
  }
  mcmc$eps_Tc[idx, country] <- eps_T_prop
  mcmc$d_c[country] <- d_prop
  warning("Cannot find likelihood increase for d_c", mcmc$meta$regions$country_code[country], ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate. = TRUE)
}


mcmc.update.U <- function(country, mcmc, ...) {
  idx <- mcmc$meta$start_c[country]:(mcmc$meta$lambda_c[country]-1)
  # exclude those extreme indices in computing loglik
  raw.outliers <- mcmc$meta$indices.outliers[[as.character(country)]]
  if (!is.null(mcmc$meta$ar.phase2) && mcmc$meta$ar.phase2) 
    raw.outliers <- sort(unique(c(raw.outliers, raw.outliers+1)))
  
  idx <- setdiff(idx, raw.outliers)
  
  z  <- (log_cond_U(mcmc$eps_Tc[idx, country],
                    mcmc$sd_Tc[idx,country],
                    mcmc$mean_eps_Tc[idx,country])
         - rexp(1))
  
  # find U_prop with g(theta_i*)>z
  # first define interval in which to search
  v <- runif(1)
  interval <- c(max(mcmc$U_c[country] - v*mcmc$meta$U.width, mcmc$meta$U.c.low[country]), 
                min(mcmc$U_c[country] + (1-v)*mcmc$meta$U.width,mcmc$meta$U.up))
  
  theta_current <-  c((mcmc$U_c[country]-mcmc$Triangle_c4[country])*exp(mcmc$gamma_ci[country,])/
                        sum(exp(mcmc$gamma_ci[country,])), mcmc$Triangle_c4[country], mcmc$d_c[country])
  theta_prop <- theta_current
  #while(TRUE) {
  for(i in 1:50) {
    U_prop <- runif(1,interval[1], interval[2])
    # keep proportions the same, just update the deltas
    theta_prop[1:3] <- theta_current[1:3]/(mcmc$U_c[country] - 
                                             mcmc$Triangle_c4[country])*(U_prop - mcmc$Triangle_c4[country])
    eps_T_prop <- get.eps.T(theta_prop, country, mcmc$meta, ...)
    if ((like <- log_cond_U(eps_T_prop, mcmc$sd_Tc[idx,country],
                   mcmc$mean_eps_Tc[idx,country])) >= z) {
      mcmc$eps_Tc[idx, country] <- eps_T_prop
      mcmc$U_c[country] <- U_prop
      return()
    } else {
      # shrink interval
      if (U_prop < mcmc$U_c[country]){
        interval[1] <- U_prop
      } else {
        interval[2] <- U_prop
      }
      if(abs(interval[1]-interval[2]) < 1e-10) break
    } # end else, interval was still ok
  }
  mcmc$eps_Tc[idx, country] <- eps_T_prop
  mcmc$U_c[country] <- U_prop
  warning("Cannot find likelihood increase for U_c", mcmc$meta$regions$country_code[country], ".\n Final interval: [", interval[1], ',', interval[2], ']\n',
          "New likelihood: ", like, ", original likelihood: ", z, immediate. = TRUE)
}