Issue using binary indicator variables as x_r in ODE system

Modeling
1

Hello everyone,

I have what ought to be a simple question. I’m trying to include multiple binary indicator variables (dummy variables) within an ODE system: B1 and B2 either take the value 0 or 1, enabling two different estimates for a parameter being fit to two groups of data; e.g., (a1*B1 + a2*B2).
I include B1 and B2 as real data x_r in the functions {} block, as well as the data {} and transformed data {} blocks (Stan code below).
The script compiles, but something is clearly going awry:

Chain 1: Rejecting initial value:
Chain 1: Error evaluating the log probability at the initial value.
Chain 1: Exception: integrate_ode_rk45: continuous data[1] is nan, but must be finite!  (in . . . at line 56)
Chain 1: Initialization between (-2, 2) failed after 100 attempts. 
Chain 1:  Try specifying initial values, reducing ranges of constrained values, or reparameterizing the model.
[1] "Error in sampler$call_sampler(args_list[[i]]) : Initialization failed."

I’ve also tried only including a single indicator variable B1 that equals a single value of 1, such that part of my ODE is essentially: a1 * 1, but this produces the same error. So I’m clearly doing something wrong in enabling x_r, but I’ve tried all possible different ways of introducing and defining x_r that I’m aware of.

Below is a simplified version of the model that produces the same error: I first simulate data in R, and then I display the Stan model. Note that in this simplified example, I don’t actually expect the parameter a1 and a2 to differ at all . . . this is just a simple example being used to diagnose the issue.

Finally, in case it’s at all useful, I also paste the full, actual Stan model I’m running, since the ODE structure and manner of calculation is a bit different (but again, even the most simplified ODE kicks up the same error).

Thank you – I very much appreciate the assistance!!

Simplified R code to simulate ODE system:

library(rstan)
library(ggplot2)
library(reshape2)

## Simulate ODE ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## params 
params <- c(a=3, 
            h=.013)

## create a list of initial conditions
R0 <- seq(50, 300, length.out=10)

## specify ODE function 
resourceLoss <- function(t, R0, params) {
  with(as.list(c(R0, params)),{
    dR <- -((a * R0)/(1 + h * a * R0))
    list(c(dR))
  })
}

## set t 
length.tlist <- 2
t.list <- seq(1, 2, length.out=length.tlist)

## run ODE 
out <- ode(y=c(R0), times=t.list, func=resourceLoss, parms=params)

## save as data frame
loss <- as.data.frame(out)

## save list without Times for later use 
dRdt.list <- round(out[,-1],3)

## create and add prefix for col names
loss.dat <- as.data.frame(dRdt.list)
colnames(loss.dat) <- paste("T", colnames(loss.dat), sep="")


## plot all dRdt 
all_dRdt <- melt(loss, id.vars="time")

tplot <- ggplot(all_dRdt, aes(time, value, col=variable,group=variable)) +
  geom_point() + geom_line(color="black") +
  xlab("Time") + ylab("Resource density") + ggtitle("Kelp remaining = R0 - dRdT") 

print(tplot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## add noise ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## reparameterized Gamma
my.rgamma<-function(n, mu, sigma){
  rgamma(n, shape = (mu^2)/(sigma^2), scale = (sigma^2)/mu)}

## set dims for matrix and set sigma 
rws=length(t.list)  ## number of rows ie data points to simulate
cls=ncol(dRdt.list)  ## number of columns ie separate FR runs 
sig=9.50

## apply rgamma to deterministic data 
mult_ODE <- function(n, dat, noise){
  noise.out <- my.rgamma(n=length(dRdt.list), mu=dat, sigma=noise)
  return(matrix(noise.out, nrow=rws, ncol=cls, byrow=FALSE))
}

noisy.dat <- round(mult_ODE(n, dat=dRdt.list, noise=sig),digits=3)

## replace noisy row 1 with deterministic R0 
noisy.dat[1,] <- dRdt.list[1,]

## create data frame for plotting, add time as a column 
noise.dat.df <- as.data.frame(noisy.dat)
noise.dat.df <- cbind(loss$time, noise.dat.df)
names(noise.dat.df)[1]<-"time"

## plot 
melt_noise <- melt(noise.dat.df, id.vars="time")

noise_plot <- ggplot(melt_noise, aes(time, value, col=variable, group=variable)) +
  geom_point() + geom_line(data=all_dRdt, aes(time, value, group=variable), col="black") +
  xlab("Time") + ylab("Resource density") + ggtitle("Resources remaining = R0 - dRdT + noise")

print(noise_plot)
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~





## configure data list for STAN ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
## extract deterministic list of initial conditions 
R0.list <- as.numeric(loss.dat[1,])

## number of columns ie ODE runs 
dim.cols <- ncol(noisy.dat)

## create binary indicator variables
ind.length = 5
B1 <- c(rep(1, length.out=ind.length), rep(0, length.out=ind.length))
B2 <- c(rep(0, length.out=ind.length), rep(1, length.out=ind.length))

## list for Stan
resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
                         t0 = t.list[1],
                         th = array(t.list[-1], dim=c(1)),
                         nLev = ncol(noisy.dat),
                         B1 = array(B1, dim=c(ncol(noisy.dat))),
                         B2 = array(B1, dim=c(ncol(noisy.dat))),
                         R0_obs = R0.list,
                         R_obs = array(noisy.dat[-1,], dim=c(1,dim.cols)))
## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

and here’s a simplified Stan model:

functions {
   real[] resourceLoss(
       real t,              // time
       real[] R,            // state 
       real[] theta,        // params
       real[] x_r, 
       int[] x_i) {
		real dRdt[x_i[1]];
		real a1 = theta[1];
		real a2 = theta[2]; 
		real h = theta[3];
		real B1 = x_r[1];   // indicator variable 1
		real B2 = x_r[2];   // indicator variable 2
		

		for (i in 1:x_i[1]) { 
	   		dRdt[i] = -(((a1 * B1) + (a2 * B2)) * R[i])/(1 + ((a1 * B1) + (a2 * B2)) * h * R[i]);
				    }
				    return dRdt;
	 	   }
}
data {
   int <lower=1> Tm;                // Tm = Tmax = T-1 (T-initial pt) = 10
   int <lower=1> nLev;              // number of diff R0 i.e. # ODEs = 
   real R0_obs[nLev];               // Intial conditions R0 list 
   real <lower=1> t0;               // starting time t = 0 
   real <lower=t0> th[Tm];          // array of length Tm i.e. 1
   real R_obs[Tm, nLev];            // simulated drift remaining
   real B1[nLev];
   real B2[nLev]; 
}
transformed data{
  real x_r[2];      
  int x_i[1];
  //real x_r[1] = B1;
  //real x_r[2] = B2;
  x_i[1] = nLev;      
}
parameters {
   real <lower=0> a1;
   real <lower=0> a2;          
   real <lower=0> h;          
   real <lower=0> sigma;      
}
transformed parameters {
   real theta[3];         
   theta[1] = a1;
   theta[2] = a2; 
   theta[3] = h;
}
model { 
   real R[Tm, nLev]; 
   real alpha[Tm, nLev];
   real beta[Tm, nLev];
 
   R = integrate_ode_rk45(resourceLoss, R0_obs, t0, th, theta, x_r, x_i);

   for (t in 1:Tm) {
     for (i in 1:nLev) { 
       alpha[t,i] = pow (R[t,i], 2) / pow (sigma, 2);
                       }
		   }

   for (t in 1:Tm) {
     for (i in 1:nLev) {
       beta[t,i] = (1 / (pow (sigma, 2) / R[t,i]));
		       }
		   }

  a1 ~ exponential(0.1);
  a2 ~ exponential(0.1);                
  h ~ exponential(0.1);                 
  sigma ~ exponential(0.1);

   for (t in 1:Tm) {
     for (i in 1:nLev) {
       R_obs[t,i] ~ gamma(alpha[t,i], beta[t,i]);
                       }
                   }  
}

In the interest of completeness, here’s the “real” system I’m wanting to fit in Stan.

dSdt = - a*S[t]*\frac{v - F[t]}{v}*q
dAdt = - a*A[t]*\frac{v - F[t]}{v}*(1 - q)
dFdt = (a*S[t]*\frac{v - F[t]}{v}*q) + (a* A[t]*\frac{v - F[t]}{v}*(1 - q)) - \frac{a*F[t]}{1 + a*b*F[t]}

And the full Stan model is below, though note I haven’t incorporated the third variable dAdt into the Stan model. But otherwise, this is the actual model I’m trying to run (I want to figure out the indicator variable before adding additional complexity).

functions {
   real[] resourceLoss(
       real t,              // time
       real[] Y,            // state 2           
       real[] theta,        // params
       real[] x_r, 
       int[] x_i) {
		real dR_dt;
		real dF_dt;
		real R = Y[1];
		real F = Y[2];
		real a = theta[1]; 
		real v1 = theta[2];
		real v2 = theta[3];
		real v3 = theta[4];
		real B1 = x_r[1];
		real B2 = x_r[2];
 		real B3 = x_r[3];
			
                dR_dt = (-(a * R) * (((v1*B1 + v2*B2 + v3*B3) - F) / (v1*B1 + v2*B2 + v3*B3)));
		dF_dt = ((a * R) * (((v1*B1 + v2*B2 + v3*B3) - F) / (v1*B1 + v2*B2 + v3*B3)));    
                return {dR_dt, dF_dt};
		    }
}
data {
   int <lower=1> Tm;                // Tm = Tmax
   int <lower=1> nLev;              // number of diff R0 i.e. # ODEs = 35
   real Y0_obs[nLev, 2];            // Intial conditions for R and F 
   real <lower=1> t0;               // starting time t = 0 
   real <lower=t0> th[Tm];          // array of length Tm 
   real R_obs[Tm, nLev];            // experimental observations
   real B1[nLev];		    // indicator variable 1: 48hrs
   real B2[nLev];		    // indicator variable 2: 96hrs
   real B3[nLev];  		    // indicator variable 3: 144hrs
}
transformed data{
  real x_r[3];
  int x_i[0];
  //real x_r[1] = B1;
  //real x_r[2] = B2;
  //real x_r[3] = B3;                
}
parameters {
   real <lower=0> a;          
   real <lower=0> v1;
   real <lower=0> v2;
   real <lower=0> v3;     
   real <lower=0> sigma;      
}
transformed parameters {
   real theta[4];
   theta[1] = a; 
   theta[2] = v1;
   theta[3] = v2;
   theta[4] = v3;
}
model { 
   real R[Tm, nLev, 2]; 
   real alpha[Tm, nLev];
   real beta[Tm, nLev];
 
   for (n in 1:nLev) {
   R[ , n, ] = integrate_ode_rk45(resourceLoss, Y0_obs[n], t0, th, theta, x_r, x_i);
   		      }

   for (t in 1:Tm) {
     for (n in 1:nLev) { 
       alpha[t,n] = pow (R[t,n,1], 2) / pow (sigma, 2);
                       }
		   }

   for (t in 1:Tm) {
     for (n in 1:nLev) {
       beta[t,n] = (1 / (pow (sigma, 2) / R[t,n,1]));
		       }
		   }
 
  a ~ exponential(0.1);                
  v1 ~ normal(50, 15); 
  v2 ~ normal(40, 25); 
  v3 ~ normal(30, 25);           
  sigma ~ exponential(0.1);

   for (t in 1:Tm) {
     for (n in 1:nLev) {
       R_obs[t,n] ~ gamma(alpha[t,n], beta[t,n]);
                       }
                   }  
}

Thanks again for any insight you all may have!!

2

It looks like you’re not actually initializing x_r?

Edit:

This does not do anything:

I guess you want something like

x_r[1] = B1;
x_r[2] = B2;

But this shouldn’t work as B1 and B2 are arrays. I guess you have to make them into reals?

3

Thanks @Funko_Unko. Yeah I commented out:

  real x_r[3];
  //x_r[1] = B1;
  //x_r[2] = B2;
  //x_r[3] = B3;

as that enabled the script to compile (but with the error described above).

I originally thought I’d indeed need to define x_r, and this is the configuration I most expected to work:

  real x_r[3];
  real x_r[1] = B1;
  real x_r[2] = B2;
  real x_r[3] = B3;

which produces:

Duplicate declaration of variable, name=x_r; attempt to redeclare as real[ ] in transformed data; previously declared as real[ ] in transformed data

If I remove real x_r[3] from the above, I get a similar error as above (it seems to view B1 and B2 as the same?:

Duplicate declaration of variable, name=x_r; attempt to redeclare as real[ ] in transformed data; previously declared as real[ ] in transformed data
 error in 'model39c43fc112c3_Fullness_IndVar' at line 39, column 18
  -------------------------------------------------
    37:   //real x_r[3];
    38:   real x_r[1] = B1;
    39:   real x_r[2] = B2;
                         ^
    40:   real x_r[3] = B3;

Note that I do first introduce real B1 = x_r[1], etc., in the functions {} block (just as I do with the thetas . . . however I tried removing these initial definitions, and instead of B1 in dR_dt, I used x_r[1], etc. . . but this produced the same issue. It seems the core issue is that Stan thinks I’m redefining previously declared variables?

I’ve tried running variations on:

 real x_r[3];
 x_r[1] = B1;
 x_r[2] = B2;
 x_r[3] = B3;

and

 //real x_r[3];
 x_r[1] = B1;
 x_r[2] = B2;
 x_r[3] = B3;

but that produces more predictable errors re: the parser expecting real, int, vector, etc.

Note that in the R script I tried explicitly defining B1, etc. as an array:

resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
                         t0 = time[1],
                         th = array(time[-1], dim=c(1)),
                         nLev = ncol(loss.dat),
                         B1 = array(B1, dim=c(dim.cols)),
                         B2 = array(B2, dim=c(dim.cols)),
                         B3 = array(B3, dim=c(dim.cols)),
                         Y0_obs = array(Y0, dim=c(dim.cols,2)),
                         R_obs = array(loss.dat[-1,], dim=c(1,dim.cols)))

and I’ve also left it “as is” (which should also be an array?):

resourceLoss.dat <- list(Tm = nrow(loss.dat)-1,
                         t0 = time[1],
                         th = array(time[-1], dim=c(1)),
                         nLev = ncol(loss.dat),
                         B1 = B1,
                         B2 = B2,
                         B3 = B3,
                         Y0_obs = array(Y0, dim=c(dim.cols,2)),
                         R_obs = array(loss.dat[-1,], dim=c(1,dim.cols)))
4

Hm, I’m a bit short on time, so I’ll just provide links to the manual.

In short, there appears to be some confusion about the Scope (computer science) - Wikipedia of variables in Stan.

Have a look here Stan Reference Manual, here Stan Reference Manual and here 7.9 Statement blocks and local variable declarations | Stan Reference Manual

Edit: I guess discourse could do a better job of naming the links there.

Edit with slightly more time:

This does not work because the first line declares x_r to be a real-valued array of size 3, while the second line attempts to declare x_r to be an array of size one and assigning to it the value of B1. Because you attempt to “declare” x_r twice in the same scope, the compiler complains about a Duplicate declaration.

If you want to assign something to the first element of x_r you have to do x_r[1] = something; The right hand side however has to be a real number, which B1 is not, at least not in the scopre you are currently in. In the transformed data block, the B1 that is seen is the one you declared in the data block, which is a real-valued array of size nLev.

None of the variables you declare in your function (block) are seen outside of the function definition, and variables from e.g. the data block are not seen inside your function definitions. Which is why you have to pass everything through x_r or x_i.

I presume your first error msg

Chain 1: Exception: integrate_ode_rk45: continuous data[1] is nan, but must be finite! (in . . . at line 56)

comes up because you did not assign any value to x_r.

1 Like
5

ahhhh okay, thank you, @Funko_Unko! That makes complete sense.

Based on your response, it sounds like the core issue is that I’m trying to assign, e.g., x_r[1] = B1, where B1 is a real-valued array of size nLev, but x_r expects a real number instead.

I want the dummy variable B1 to be of length nLev because nLev is the number of separate runs of the ODE system (each run w/ different initial conditions, and where the data for initial conditions is real R0_obs[nLev] in the data block). Thus, I want each ODE run to be associated with a value of B1

I’m not precisely sure how to do this, though perhaps I can modify the loop that I already have in place for integrate_ode_rk45(). This is from the second model in my first post, and here I’ve added the [n] to x_r[n]). Tm is the array of time points we want solutions for, and the 2 is for a second state variable (but after thinking about this, I’m now somewhat doubtful this x_r[n] would work if I have multiple indicator variables, i.e. B1, B2, B3:

real R[Tm, nLev, 2];
   for (n in 1:nLev) {
   R[ , n, ] = integrate_ode_rk45(resourceLoss, Y0_obs[n], t0, th, theta, x_r[n], x_i);
   		      }

In either case, I still need to pass a real-valued array of length nLev as x_r . . . or find some other way of supplying that information. Are there workarounds to supply an array to x_r, or perhaps is there another way of providing that information within integrate_ode_rk45()?

Thank you, I really appreciate the help!!

6

I think something like this should be possible:

real R[Tm, nLev, 2];
   for (n in 1:nLev) {
   R[ , n, ] = integrate_ode_rk45(resourceLoss, Y0_obs[n], t0, th, theta, {B1[n],B2[n}, x_i);
   		      }
1 Like
7

Thank you, @Funko_Unko, I really appreciate your help! That did it!

Instead of attempting to explicitly define any of the B1, B2, etc., in the transformed parameter block, I now simply state real x_r[3] (and where B1, B2, & B3 are still defined as real B1[nLev] in the data block). Then, per your suggestion, the indicator variables are called within the integrate_ode_rk45{} loop:

for (n in 1:nLev) {
R[ , n, ] = integrate_ode_rk45(resourceLoss, Y0_obs[n], t0, th, theta, {B1[n], B2[n], B3[n]}, x_i);
   		                }

The model is running now and will take a little while, but it has advanced past the stage where all errors pertaining to this issue had manifested previously. So, I think I’m all set :-)

Again, thank you!! I’m sure these are fairly obvious questions, and I’m simply new to Stan (and programming in general), and thus I still get caught on what are surely simple issues.

Cheers!
Zach

8

I’m glad that it works now.

Don’t forget to check for other problems: Brief Guide to Stan’s Warnings