Cant run Stan Code

Modeling
1

Please share your Stan program and accompanying data if possible.

When including Stan code in your post it really helps if you make it as readable as possible by using Stan code chunks (```stan) with clear spacing and indentation. For example, use

model {
  vector[N] mu = alpha + beta * x;
  y ~ normal(mu, sigma); 
}

instead of

model{
vector[N] mu = alpha+beta*x;
y~normal(mu,sigma);
}

To include mathematical notation in your post put LaTeX syntax between two $ symbols, e.g.,
p(\theta | y) \propto p(\theta) p(y | \theta).

Don’t forget to add relevant tags to your topic (top right of this form) for application area and/or class of models you work with.

Hello i am trying to fit my stan model but it is giving me the following error

Error in mod$fit_ptr() : 
  Exception: variable does not exist; processing stage=data initialization; variable name=infection_times; base type=double  (in 'model3c29ba681a_DSA' at line 25)

failed to create the sampler; sampling not done

Stan model 'DSA' does not contain samples.
NULL
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.
Stan model 'DSA' does not contain samples.

Error in hist.default(beta_sample): 'x' must be numeric
Traceback:

1. hist(beta_sample)

2. hist.default(beta_sample)

3. stop("'x' must be numeric")

Below is my stan code

functions{
  real[] SIR(real t, real[] y, real[] parms, real[] x_r, int[] x_i){
    real beta = parms[1];
    real gamma = parms[2];
    real dydt[3];
    dydt[1] = -beta * y[1] * y[2];
    dydt[2] = beta * y[1] * y[2] - gamma * y[2];
    dydt[3] = gamma * y[2];
    return dydt;
  }
  
  real[] SIR_diffusion(real t, real[] y, real[] parms, real[] x_r, int[] x_i){
    real beta = parms[1];
    real gamma = parms[2];
    real dydt[3];
    dydt[1] = -beta * y[1] * y[2];
    dydt[2] = beta * y[1] * y[2];
    dydt[3] = 0;
    return dydt;
  }
}

data{
  int<lower=0> N; 
  real<lower=0.0> infection_times[N];
}


transformed data {
    real x_r[0];
    int x_i[0];
}

parameters {
    real<lower=0.0> beta;
    real<lower=0.0> gamma;
    real<lower=0, upper=1.0> rho;
}

transformed parameters{
  real<lower=0.0> R0 = beta/gamma;
}

model{
  real parms[3];
  real ic[3];
  real sol[N,3];
  real t0;
  real Smax;
  real tau_T;
  
  parms[1] = beta;
  parms[2] = gamma;
  parms[3] = rho;
  
  ic[1] = 1.0;
  ic[2] = rho;
  ic[3] = 0.0;
  
  t0 = 0.0; 
  
  sol = integrate_ode_rk45(SIR,ic,t0,infection_times,parms,x_r,x_i);
  Smax = sol[N,1];
  tau_T = 1.0 - Smax;
  
  for (i in 1:N){
    target += log((beta * sol[i,1] * sol[i,2])/tau_T);
  }
  
  target += gamma_lpdf(beta | 2, 2) + gamma_lpdf(gamma | 2, 2) + beta_lpdf(rho | 1, 1);
  }

and here is my R code

#options(warn=-1, message =-1)
library(rstan); 
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)
library(ggplot2); 



#### Import data 
sellke.data <- read.csv("infection_times.csv", header =  TRUE)
infection_times <- sellke.data[sellke.data[,1] > 0, 1,drop=FALSE]$infection_times


## Stan configurations
N <-1000 #sample size
nIter <-10000 #number of iterations of Hamiltonian Monte Carlo 

## Draw random samples of infection times

if (length(infection_times) < N){
  infection_times <- sort(unique(infection_times), decreasing = FALSE)
  N <- length(infection_times)
} else {
  infection_times <- sample(unique(infection_times), size = N, replace = FALSE)
  infection_times <- sort(unique(infection_times), decreasing = FALSE)
  N <- length(infection_times)
}

stan_data <- list(N=N, beta=beta, gamma=gamma, infection_times = infection_times)

fit <- stan(
  file = "DSA.stan",
  data = stan_data,
  chains = 1,
  iter = nIter,
 refresh = 20000
)

stan_summary_fit<-summary(fit)
print(stan_summary_fit)


#Trace plots
traceplot(fit, "beta")
traceplot(fit, "gamma")
traceplot(fit, "rho")

# Posterior histograms
beta_sample <- extract(fit)$beta
gamma_sample <- extract(fit)$gamma 
rho_sample <- extract(fit)$rho
R0_sample <- extract(fit)$R0



hist(beta_sample)
hist(gamma_sample)
hist(rho_sample)
hist(R0_sample)
2

I did some edits to your post to make it easier to read. A couple other things that will help folks help you out.

  1. The OS and version of Stan you are using.
  2. Snippet of your data or some simulated data

Thanks! And welcome.