Hi
I’m trying to understand why with the following model cumulative_R has a much narrower posterior than cumulative_ncp, given that I would expect them to be identical.
In this model there are two rates, things that vary at rate N and things that vary at rate day;
cumulative_ncp is the sum of the variable delta varying at N.
while cumulative_R is the sum of delta (as in this cut-down version of the model delta = R) varying at day. both are then scaled to the same units.
In my test case where the sum of delta x 86400 is known to be (-)700.
The final value of cumulative_R is 697.7598-701.9972 (5-95% quantile range) while
cumulative_ncp is -689.3504-704.6502.
data {
int<lower=0> N;
int<lower=0> N_day;
int<lower=0> day[N];
real C[N];
real kw[N];
real B[N];
real hm[N];
real dhdt[N];
real Csat[N];
real dt;
}
parameters {
real<lower=0> sigma;
real<lower=0> Cb[N_day];
real<lower=0> R[N_day];
}
transformed parameters {
real C_hat[N-1];
real delta[N-1];
for(n in 1:N-1){
delta[n] = 0 - R[day[n]];
C_hat[n] = (exp(-(dt*(kw[n] + dhdt[n]))/hm[n]) * (dhdt[n] * (Cb[day[n]] * (exp((dt*(kw[n] + dhdt[n]))/hm[n]) - 1) + C[n]) + (delta[n]) * (exp((dt * (kw[n] + dhdt[n]))/hm[n]) - 1) + kw[n] * ( (Csat[n] * ((B[n]) + 1)) * (exp((dt * (kw[n] + dhdt[n]))/hm[n]) - 1) + C[n])))/(kw[n] + dhdt[n]);
}
}
model {
sigma ~ normal(0, 1);
R ~ normal(0, 1);
Cb ~ normal(200, 10);
for(n in 2:N){
C[n] ~ normal(C_hat[n-1], sigma);
}
}
generated quantities{
vector[N-1] cumulative_ncp;
vector[N_day] cumulative_R;
for(n in 1:N-1){
cumulative_ncp[n] = (delta[day[n]]); // mmol m-2 s-1
}
cumulative_ncp = cumulative_sum(ncp * dt); // mmol m-2 day-1
for(d in 1:N_day){
cumulative_R[d] = (R[d]) * 86400; // mmol m-2 day-1
}
cumulative_R = cumulative_sum(cumulative_R);
}