data{
    int count[128];
    vector[128] time;
}
parameters{
    real logK1_std;
    real<lower=0> r1_std;
    real<lower=0, upper = 1> L01_prop;
    real logK2_std;
    real<lower=0> neg_r2_std;
    real<lower=0, upper = 1> L02_prop;
    real<lower=0> theta_std;
}
transformed parameters{
    real<lower = 0> K1 = exp(log(10) + logK1_std * 3 * log(10));
    real<lower = 0> K2 = exp(log(10) + logK2_std * 3 * log(10));
    real<lower = 0> L01 = K1 * L01_prop;
    real<lower = 0> L02 = K2 * L02_prop;
    real<lower = 0> r1 = r1_std / 1;
    real<upper = 0> r2 = -neg_r2_std / 1;
    real<lower=0> theta = theta_std / 0.5;
}
model{
    vector[128] lambda;
    vector[128] alpha;
    vector[128] beta;
    theta_std ~ exponential( 1 );

    //growth model
    r1_std ~ exponential( 1 );
    L01_prop ~ beta(1, 100);
    logK1_std ~ std_normal();

    //decay model
    neg_r2_std ~ exponential( 1 );
    L02_prop ~ beta(100, 1);
    logK2_std ~ std_normal();

    //mathemagically combine them
    for ( i in 1:128 ) {
        lambda[i] = K1 / (1 + ((K1 - L01) / L01) * exp(-r1 * time[i])) +
                    K2 / (1 + ((K2 - L02) / L02) * exp(-r2 * time[i]));
        alpha[i] = lambda[i] / theta;
        beta[i] = 1 / theta;
    }

    //evaluate likelihood
    count ~ neg_binomial( alpha , beta );
}
generated quantities{
	//none for now
}