Strange new error when running code that has always worked before

Modeling
1

I have just tried to rerun in Rstudio some code that always worked previously but I get a strange error message

Error in cpp_object_initializer(.self, .refClassDef, ...) : 
  could not find function "cpp_object_initializer"

That error message was also found herehttps://discourse.mc-stan.org/t/gaussian-process-es-gp-exp-quad-cov-not-found/6481/7 but I don’t see a solution there.

I am using: rstan 2.18.2, StanHeaders 2.18.0 in Rstudio 1.1.453 with R 3.5.1

Here is my Stan code, but it has worked multiple( = thousands of) times without any problem).

The purpose of the code is to generate a variance-covariance matrix using a somewhat complex kernel function.

functions {


  real gammaseq(real theta,
                    real phi,
                      real alpha,
                      real beta,
                      real sigma_sq,
                      real b){
            real term1 =   log(sigma_sq) +
                          2 * (alpha * log(beta) - lgamma(alpha));


            real term2 =  (alpha - 1) * (log(theta) + log(phi));

            real term3 =  -beta * (theta + phi);

            real term4 =  -inv(2 * b) * square(theta - phi);

            return( exp(term1 + term2 + term3 + term4));
                      }


    // It is a function of x
  real first_test_integral( real x,
                      real ymax,
                      int polyorder,
                      real[] roots,
                      real[] weights,
                      real alpha,
                      real beta,
                      real sigma_sq,
                      real b){

        real cum[polyorder]; // accumulator for result

        for ( i in 1:polyorder){
          cum[i] =   weights[i] * gammaseq(
            x,
            0.5 * ymax * roots[i] + 0.5 * ymax,
                        alpha, beta, sigma_sq, b);
          }

              return(0.5 * ymax * sum(cum)) ;
   }
  real second_test_integral(real xmax,
                      real ymax,
                      int polyorder,
                      real[] roots,
                      real[] weights ,
                      real alpha,
                      real beta,
                      real sigma_sq,
                      real b){

        real cum[polyorder]; // accumulator for result

    for ( i in 1:polyorder){

      cum[i] = weights[i] * first_test_integral(
        0.5 * xmax * roots[i] + 0.5 * xmax,
        ymax,
        polyorder,
        roots,
        weights,
                        alpha, beta, sigma_sq, b);
          }

        return(0.5 * xmax * sum(cum));
        }


}
data {
  int N;
  int M;
  int polyorder;
  vector[M] x ;
  vector[N] y ;
  real<lower = 0> alpha;
  real<lower = 0> beta;
  real<lower = 0> sigma_sq;
  real<lower = 0> b;
  real roots[polyorder];
  real weights[polyorder];



}

model {}
generated quantities{
  matrix[M, N] K;
  for(j in 1:N){
    for(i in 1:M){
      K[i, j] = sigma_sq * second_test_integral(x[i],
                               y[j],
                               polyorder,
                               roots,
                               weights,
                        alpha, beta, sigma_sq, b);
    }
  }

}

The R code that drives it is this function:

kxx_gamma <- function(X,
                      Y,
                      alpha,
                      beta,
                      sigma_sq,
                      b,
                      roots,
                      weights,
                      stanmodl ){
  M <- length(X)
  N <- length(Y)
  standata <- list(x = X,
                   y = Y,
                   polyorder = length(roots),
                   M = M,
                   N = N,
                   sigma_sq = sigma_sq,
                   alpha = alpha,
                   beta = beta,
                   b = b,
                   roots = roots,
                   weights =  weights)
  fit <-purrr::quietly(rstan::sampling)( stanmodl,
                         data = standata,
                         chains = 1,
                         iter = 1,
                         algorithm = "Fixed_param"
  )$result
  K <- rstan::extract(fit, pars = "lp__", include = FALSE)
  KK <- matrix(unlist(K), ncol = N, nrow = M)
  return(KK)
}

The error arises instantly on the sampling call.

2

rstan::sampling does not work for longstanding but never fully understood reasons. You have to do

library(rstan)
sampling(...)
3

OK. Many thanks for super-fast reply.

Doesn’t that cause a problem if I wanted to include that code in an R package?

4

Inside another package it works if you depend on Rcpp and import sampling from rstan.