Smoothed logistic growth with time-varying carrying capacity.

corr_mod.stan

data {
  int<lower=0> num_obs; // number of observations
  real y[num_obs]; // Observed value
  int time[num_obs]; // time
  
}
parameters {
  //weakly uninformative priors, K on log base 10 scale
  vector[100] log_K; 
  real<lower=0,upper=100> p0; 
  real r; 
  real<lower=0,upper=5>  sigma;
}

transformed parameters {
  vector[100] K; 
  K=exp(log_K);
}

model {
  
  real y_mean[num_obs] ; 
  vector[num_obs] y_sigma;
  
  log_K[1] ~ normal(3, .2);
  for( i in 2:100){
    log_K[i] ~ normal(log_K[i-1], .1);
  }
  
  
  for(i in 1:num_obs){
    y_mean[i] = K[time[i]]*1/(1+((K[time[i]]-p0)/p0)*exp(-r*time[i]) ) ;
  }
  
  for(i in 1:num_obs){
    y_sigma[i] = sigma*y_mean[i] ;
  }
  
  y ~normal(y_mean, y_sigma) ;
}

run_stan.R

library(rstan)

set.seed(717)

##Generate data
##use functional form from ecology: K*1/(1+((K-p0)/p0)*exp(-r*t) )
##see https://en.wikipedia.org/wiki/Logistic_function#In_ecology:_modeling_population_growth
##K - long term maximum of S curve - settling to 100
## p0 initial value - setting to 1
## r - rate of growth --setting to .1

K=100
p0=1
r=.1

ti=1:100
z=K*1/(1+((K-p0)/p0)*exp(-1*r*ti) )

##sample z with heteroskedasctic 20% error, note in STAN we will estimate this parameter
y=rnorm(n=length(z), mean=z, sd=z*.2)

## only use first 40 time points
stan_dat <- list(
  num_obs = 40,
  y = y[1:40],
  time=ti[1:40]
)

#### --------    Run time-varying K (all correlated)  -------      ####
corr_mod = stan_model("corr_mod.stan")

corr_fit <- sampling(corr_mod, 
                     data = stan_dat,
                     iter = 5000, thin=1, chains = 1)

corr_res=extract(corr_fit)

#### --------    Run constant-K model   -------      ####
van_mod = stan_model("vanilla_mod.stan")


van_fit <- sampling(van_mod, 
                     data = stan_dat,
                     iter = 5000, thin=1, chains = 1)

van_res=extract(van_fit)

## projection
par(mfrow=c(1,2))

plot(ti,z, ylim=c(0,200), type='l')
for(i in 1:1000){
  lines(ti, 
        corr_res$K[i,]*1/(1+((corr_res$K[i,]-corr_res$p0[i])/corr_res$p0[i])*exp(-1*corr_res$r[i]*ti) ),
        col='lightgray')
}
lines(ti, z, col='black')
points(ti, y, pch=20, col='black')
abline(v=40, lty=2)


plot(ti,z, ylim=c(0,200), type='l')
for(i in 1:1000){
  lines(ti, 
        van_res$K[i]*1/(1+((van_res$K[i]-van_res$p0[i])/van_res$p0[i])*exp(-1*van_res$r[i]*ti) ),
        col='lightgray')
}
lines(ti, z, col='black')
points(ti, y, pch=20, col='black')
abline(v=40, lty=2)

vanilla_mod.stan

data {
  int<lower=0> num_obs; // number of observations
  real y[num_obs]; // Observed value
  real time[num_obs]; // time
}
parameters {
//weakly uninformative priors, K on log base 10 scale
  real<lower=0,upper=4> log_K; 
  real<lower=0,upper=100> p0; 
  real<lower=-0.5,upper=0.5> r; 
  real<lower=0,upper=5>  sigma;
}

transformed parameters {
  real K; 
  K=10^log_K;
}

model {
  real y_mean[num_obs] ; 
  vector[num_obs] y_sigma;
  
  
  for(i in 1:num_obs)
  		y_mean[i] = K*1/(1+((K-p0)/p0)*exp(-r*time[i]) ) ;
      
  for(i in 1:num_obs)
  		y_sigma[i] = sigma*y_mean[i] ;
  
  y ~normal(y_mean, y_sigma) ;
}