Rcpp function within R loop crashes, whereas outside of the loop it doesn't

Question

I wrote the following code in Rcpp

  #include <RcppArmadilloExtensions/sample.h>
  #include <random>
  #include <iostream>
  #include <cmath>
  #include <vector>
  #include <numeric>
  #include <algorithm>
  #include <stdio.h>    
  #include <math.h>
  #include<Rmath.h>
  using namespace Rcpp;
// [[Rcpp::export]]
double StrausProcess_Unc(NumericVector xi, double MX, double MN, double a, double d, double b){
  if (sum(xi > MX) == 0 & sum(xi < MN) == 0) {
    int K_dot = xi.size();
    double term_2 = 0;
    if (K_dot > 1) {
      for (int i = 0; i < (K_dot - 1); ++i) {
        for (int j = (i + 1); j < K_dot; ++j) {
          term_2 += (abs(xi[i] - xi[j]) <= d) * log(a);
        }
      }
    }
    return K_dot * log(b) + term_2;
  } else {
    for (int i = 0; i < xi.size(); ++i) {
      xi[i] = std::max(std::min(xi[i], MX), MN);
    }
    int K_dot = xi.size();
    double term_2 = 0;
    if (K_dot > 1) {
      for (int i = 0; i < (K_dot - 1); ++i) {
        for (int j = (i + 1); j < K_dot; ++j) {
          term_2 += (abs(xi[i] - xi[j]) <= d) * log(a);
        }
      }
    }
    return K_dot * log(b) + term_2;
  }
}


// [[Rcpp::export]]
NumericVector Birth_Death_Sim_hard_sample_2(int niter, NumericVector mu_init, double q_move, double q_birth,
                                          double a, double d, double beta, double MX, double MN,
                                          double Upper, double Lower){
  double num;
  double den;
  int old_n;
  int n = mu_init.size();
  NumericVector mu_cand(n);
  NumericVector mu_old(n);
  mu_old = clone(mu_init);
  double m_birth;
  NumericVector m_birth_1(2); //vector to be used only when we reach length 1
  int j_star; //position for birth sample
  int j_star_1; //position of component to be killed
  for(int iter=1; iter<niter; ++iter){
    //move change current means
    old_n = mu_old.size();
    if(R::runif(0,1) <= q_move){
      for(int k=0; k<old_n; ++k){
        mu_cand[k] = R::rlnorm(log(mu_old[k]),0.1);
      }
      num = StrausProcess_Unc(mu_cand,  MX,  MN,  a,  d,  beta);

      den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta);
      if(R::runif(0,1)<= exp(num-den)){
        mu_old = clone(mu_cand);
      }else{
        mu_old = clone(mu_old);
      }
    }else{
      if(old_n==1){// give birth
        m_birth = R::runif(Lower,Upper);//R::rlnorm(mu_0,sigma_0);
        m_birth_1[0] = mu_old[0];
        m_birth_1[1] = m_birth;
        num = StrausProcess_Unc(m_birth_1,  MX,  MN,  a,  d,  beta);
        den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta)  + R::dunif(m_birth,Lower,Upper,1);
        if(R::runif(0,1) <= exp(num-den)){
          mu_old = clone(m_birth_1);
        }else{
          mu_old = clone(mu_old);
        }
      }else{
         if(R::runif(0,1) <= q_birth & old_n<100){//give birth
          j_star = RandInt(old_n+1);
          NumericVector mu_cand_i(old_n+1);
          mu_cand_i[j_star] = R::runif(Lower,Upper);//R::rlnorm(mu_0,sigma_0);
          if(j_star==0){
            Range r(j_star + 1, old_n);
            mu_cand_i[r] = clone(mu_old);
          }else if(j_star== old_n){
            Range r(0, old_n - 1);
            mu_cand_i[r] = clone(mu_old);
          }else{
             Range r(0, j_star-1);
             Range u(j_star+1,old_n);
             Range x(j_star,old_n-1);
             mu_cand_i[r] = mu_old[r];
             mu_cand_i[u] =  mu_old[x];
          }
          num = StrausProcess_Unc(mu_cand_i,  MX,  MN,  a,  d,  beta);
          den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta) + R::dunif(mu_cand_i[j_star],Lower,Upper,1);
          if(R::runif(0,1) <= exp(num-den)){
            mu_old = clone(mu_cand_i);
          }else{
            mu_old = clone(mu_old);
          }
        }else{//give death
          j_star_1 = RandInt(old_n);
          NumericVector mu_cand_j(old_n-1);
          if(j_star_1==0){
            Range r(j_star_1+1, old_n-1);
            mu_cand_j = mu_old[r];
          }else if(j_star_1 == (old_n-1) ){
            Range r(0, old_n-2);
            mu_cand_j = mu_old[r];
          }else{
             Range r(0, j_star_1-1);
             Range u(j_star_1+1,old_n-1);
             Range x(j_star_1,old_n-2);
             mu_cand_j[r] = mu_old[r];
             mu_cand_j[x] =  mu_old[u];
          }
          num = StrausProcess_Unc(mu_cand_j,  MX,  MN,  a,  d,  beta) + R::dunif(mu_old[j_star_1],Lower,Upper,1);
          den = StrausProcess_Unc(mu_old,  MX,  MN,  a,  d,  beta);
          if(R::runif(0,1)<=exp(num-den)){
            mu_old = clone(mu_cand_j);
          }else{
            mu_old = clone(mu_old);
          }
        }
      }
    }
  }
  return(mu_old);

}

However, my problem is the following, when I call the function Birth_Death_Sim_hard_sample_2 within a loop in R after a while it is always crashes.

For example, you can run

sz = c()
for(i in 1:20000){
sz[i] =   length( Birth_Death_Sim_hard_sample_2(50, c(0.5, 5.0, 10.0, 250.0, 500.0), 0.5, 0.5,
                                1e-300,  35, 5, 2400,  0,
                                242, 2)
}

At least on my machine it never manages to finish the loop. Am I doing something wrong? Because when I compile the code I get no error and also outside the loop I can call the Birth_Death_Sim_hard_sample_2 repeteadly without crashing.