Sampler fails when using custom hybrid likelihood function written with a mixture of R and Rcpp about drjacoby HOT 12 OPEN

// [[Rcpp::export]]
NumericVector likelihood_fast(const NumericVector &expected, const NumericVector &obs, const NumericVector &pars){
    int total_cts = expected.size();
    NumericVector ret(total_cts);
    const double sd = pars["error"];
    const double max_ct = 40.0;
    const double log_const = log(0.5);
    const double den = sd*M_SQRT2;
    
    for(int i = 0; i < total_cts; ++i){
        if(obs[i] < max_ct && obs[i] >= 0.0){
            ret[i] = -log(sd*sqrt(2*M_PI)) - 0.5*pow((obs[i] - expected[i])/sd, 2);
        } else if(obs[i] >= max_ct) {
            ret[i] = log_const + log(erfc((max_ct - expected[i])/den));
        } else {
            ret[i] = log_const + log(1.0 + erf((0.0 - expected[i])/den));
        }
    }
    return(ret);
} 

// [[Rcpp::export]]
double ct_likelihood_fast_tmp(double x){
    return(x);
} 

library(drjacoby)
library(Rcpp)


cppFunction("NumericVector dnorm_cpp(const NumericVector &expected, const double &obs, const NumericVector &pars){
    int n_ret = expected.size();
    NumericVector ret(n_ret);
    const double sd = pars[0];
    const double max_ct = 40.0;
    const double log_const = log(0.5);
    const double den = sd*M_SQRT2;
    
    for(int i = 0; i < n_ret; ++i){
        ret[i] = -log(sd*sqrt(2*M_PI)) - 0.5*pow((obs - expected[i])/sd, 2);
    }
    return(ret);
}")

set.seed(1)

# define true parameter values
mu_true <- 3
sigma_true <- 2

# draw example data
data_list <- list(x = rnorm(10, mean = mu_true, sd = sigma_true))

# define parameters dataframe
df_params <- define_params(name = "mu", min = -10, max = 10,
                           name = "sigma", min = 0, max = Inf)

# define log-likelihood function
r_loglike <- function(params, data, misc) {
    
    # extract parameter values
    mu <- params["mu"]
    sigma <- params["sigma"]
    
    # calculate log-probability of data
    ret <- sum(dnorm(data$x, mean = mu, sd = sigma, log = TRUE))
    
    # return
    return(ret)
}

# define log-likelihood function
r_loglike_cpp <- function(params, data, misc) {
    
    # extract parameter values
    mu <- params["mu"]
    sigma <- params["sigma"]
    
    # calculate log-probability of data
    ret <- sum(dnorm_cpp(data$x, mu, sigma))
    
    # return
    return(ret)
}


# define log-prior function
r_logprior <- function(params, misc) {
    
    # extract parameter values
    mu <- params["mu"]
    sigma <- params["sigma"]
    
    # calculate log-prior
    ret <- dunif(mu, min = -10, max = 10, log = TRUE) +
        dlnorm(sigma, meanlog = 0, sdlog = 1.0, log = TRUE)
    
    # return
    return(ret)
}

## Run example
mcmc <- run_mcmc(data = data_list,
                 df_params = df_params,
                 loglike = r_loglike,
                 logprior = r_logprior,
                 burnin = 1e3,
                 samples = 1e3,
                 pb_markdown = TRUE)
plot_par(mcmc, show="mu")

## Run example with custom Cpp dnorm
mcmc_cpp <- run_mcmc(data = data_list,
                 df_params = df_params,
                 loglike = r_loglike_cpp,
                 logprior = r_logprior,
                 burnin = 1e3,
                 samples = 1e3,
                 pb_markdown = TRUE)
plot_par(mcmc_cpp, show="mu")

library(drjacoby)

# Source an empty cpp function
Rcpp::cppFunction("void break_mcmc(){}", verbose=TRUE)

# define log-likelihood function
ll <- function(params, data, misc) {
  # Option to call the empty cpp function
  if(misc$broken){
    break_mcmc()
  }
  # calculate log-probability of data
  ret <- sum(dnorm(data$x, mean = params["mu"], sd = 2, log = TRUE))
  # return
  return(ret)
}

# define log-prior function
lp <- function(params, misc) {
  return(0)
}

# define true parameter values
mu_true <- 3
sigma_true <- 2
# draw example data
set.seed(1234)
data_list <- list(x = rnorm(10, mean = mu_true, sd = sigma_true))
# define parameters dataframe
df_params <- define_params(name = "mu", min = -10, max = 10, init = 3)

# Working
set.seed(1234)
mcmc_cpp <- run_mcmc(data = data_list,
                     df_params = df_params,
                     loglike = ll,
                     logprior = lp,
                     burnin = 500,
                     samples = 500,
                     chains = 1,
                     silent = TRUE,
                     misc = list(broken = FALSE))
# Returns MCMC chain for mu that looks fine:
plot(mcmc_cpp$output$mu, t = "l")

# Not working
set.seed(1234)
mcmc_cpp2 <- run_mcmc(data = data_list,
                     df_params = df_params,
                     loglike = ll,
                     logprior = lp,
                     burnin = 500,
                     samples = 500,
                     chains = 1,
                     silent = TRUE,
                     misc = list(broken = TRUE))
# Returns MCMC chain for mu that does not look fine (either monotonically increases or decreases):
plot(mcmc_cpp2$output$mu, t = "l")