Hello there,
I am performing a Bayesian model for analysis of estuaries metabolism. I am editing the original code by
which allows to customize the model equations in Stan. I tried updating the ER function with a power function (ER(t) = ER20 * 1.045 ^ (temp.water(t) - 20) instead of a linear function, and I got the following error.
I also tried using 1.045^ ( temp_water[i] - 20) but I also got syntax error.
Any idea how to go about this?
My R code looks like this:
mm_saturator_ER<-
mm_name(‘bayes’, GPP_fun=‘satlight’, ER_fun=‘q10temp’,pool_K600=‘none’) %>%
specs(alpha_meanlog=.4,alpha_sdlog=.2,Pmax_mu=6,Pmax_sigma=.3, day_start = 4,day_end = 28,n_cores=4, n_chains=4, burnin_steps=500, saved_steps=500) %>%
metab(dat)
######################
SYNTAX ERROR, MESSAGE(S) FROM PARSER:
No matches for:
pow(vector)
Available argument signatures for pow:
pow(real, real)
error in ‘model506852d32a8b_b_np_oipi_tr_psrqkm’ at line 76, column 52
74: for(i in 1:n) {
75: GPP_inst[i] = Pmax .* tanh(light[i] .* alpha ./ Pmax);
76: ER_inst[i] = ER20 *1.045 .*pow(temp_water[i]-20);
^
77: KO2_inst[i] = K600_daily .* KO2_conv[i];
Timing stopped at: 0.08 0 0.08
Warning message:
In metab_fun(specs = specs, data = data, data_daily = data_daily, :
Modeling failed
Errors:
failed to parse Stan model ‘b_np_oipi_tr_psrqkm’ due to the above error.
Please share your Stan program and accompanying data if possible.
When including Stan code i
n 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
// b_np_oipi_tr_psrqkm.stan
data {
// Parameters of priors on metabolism
real alpha_meanlog;
real<lower=0> alpha_sdlog;
real<lower=0> Pmax_mu;
real<lower=0> Pmax_sigma;
real<lower=0> ER_mu;
real<lower=0> ER_sigma;
real<lower=0> ER20_mu;
real<lower=0> ER20_sigma;
real K600_daily_meanlog;
real<lower=0> K600_daily_sdlog;
// Error distributions
real<lower=0> err_obs_iid_sigma_scale;
real<lower=0> err_proc_iid_sigma_scale;
// Data dimensions
int<lower=1> d; // number of dates
real<lower=0> timestep; // length of each timestep in days
int<lower=1> n24; // number of observations in first 24 hours per date
int<lower=1> n; // number of observations per date
// Daily data
vector[d] DO_obs_1;
// Data
vector[d] DO_obs[n];
vector[d] DO_sat[n];
vector[d] light[n];
vector[d] depth[n];
vector[d] temp_water[n];
vector[d] KO2_conv[n];
}
parameters {
vector[d] alpha_scaled;
vector[d] Pmax;
vector[d] ER;
real ER20;
vector<lower=0>[d] K600_daily;
real<lower=0> err_obs_iid_sigma_scaled;
real<lower=0> err_proc_iid_sigma_scaled;
vector[d] DO_mod[n];
}
transformed parameters {
real<lower=0> err_obs_iid_sigma;
vector[d] DO_mod_partial_sigma[n];
real<lower=0> err_proc_iid_sigma;
vector<lower=0>[d] alpha;
vector[d] GPP_inst[n];
vector[d] ER_inst[n];
vector[d] KO2_inst[n];
vector[d] DO_mod_partial[n];
// Rescale error distribution parameters
err_obs_iid_sigma = err_obs_iid_sigma_scale * err_obs_iid_sigma_scaled;
err_proc_iid_sigma = err_proc_iid_sigma_scale * err_proc_iid_sigma_scaled;
// Rescale select daily parameters
alpha = exp(alpha_meanlog + alpha_sdlog * alpha_scaled);
// Model DO time series
// * trapezoid version
// * observation error
// * IID process error
// * reaeration depends on DO_mod ERT is ER20
// Calculate individual process rates
for(i in 1:n) {
GPP_inst[i] = Pmax .* tanh(light[i] .* alpha ./ Pmax);
ER_inst[i] = ER20 *1.045 .*pow(temp_water[i]-20);
KO2_inst[i] = K600_daily .* KO2_conv[i];
}
// DO model
DO_mod_partial[1] = DO_obs_1;
DO_mod_partial_sigma[1] = err_proc_iid_sigma * timestep ./ depth[1];
for(i in 1:(n-1)) {
DO_mod_partial[i+1] =
DO_mod[i] .*
(2.0 - KO2_inst[i] * timestep) ./ (2.0 + KO2_inst[i+1] * timestep) + (
(GPP_inst[i] + ER_inst[i]) ./ depth[i] +
(GPP_inst[i+1] + ER_inst[i+1]) ./ depth[i+1] +
KO2_inst[i] .* DO_sat[i] + KO2_inst[i+1] .* DO_sat[i+1]
) .* (timestep ./ (2.0 + KO2_inst[i+1] * timestep));
for(j in 1:d) {
DO_mod_partial_sigma[i+1,j] = err_proc_iid_sigma *
sqrt(pow(depth[i,j], -2) + pow(depth[i+1,j], -2)) .*
(timestep / (2.0 + KO2_inst[i+1,j] * timestep));
}
}
}
model {
// Independent, identically distributed process error
for(i in 1:n) {
DO_mod[i] ~ normal(DO_mod_partial[i], DO_mod_partial_sigma[i]);
}
// SD (sigma) of the IID process errors
err_proc_iid_sigma_scaled ~ cauchy(0, 1);
// Independent, identically distributed observation error
for(i in 1:n) {
DO_obs[i] ~ normal(DO_mod[i], err_obs_iid_sigma);
}
// SD (sigma) of the observation errors
err_obs_iid_sigma_scaled ~ cauchy(0, 1);
// Daily metabolism priors
alpha_scaled ~ normal(0, 1);
Pmax ~ normal(Pmax_mu, Pmax_sigma);
ER ~ normal(ER_mu, ER_sigma);
ER20 ~ normal(ER20_mu, ER20_sigma);
K600_daily ~ lognormal(K600_daily_meanlog, K600_daily_sdlog);
}
generated quantities {
vector[d] err_obs_iid[n];
vector[d] err_proc_iid[n-1];
vector[d] GPP;
vector[d] ER;
vector[n] DO_obs_vec; // temporary
vector[n] DO_mod_vec; // temporary
vector[d] DO_R2;
for(i in 1:n) {
err_obs_iid[i] = DO_mod[i] - DO_obs[i];
}
for(i in 2:n) {
err_proc_iid[i-1] = (DO_mod_partial[i] - DO_mod[i]) .* (err_proc_iid_sigma ./ DO_mod_partial_sigma[i]);
}
for(j in 1:d) {
GPP[j] = sum(GPP_inst[1:n24,j]) / n24;
ER[j] = sum(ER_inst[1:n24,j]) / n24;
// Compute R2 for DO observations relative to the modeled, process-error-corrected state (DO_mod)
for(i in 1:n) {
DO_mod_vec[i] = DO_mod[i,j];
DO_obs_vec[i] = DO_obs[i,j];
}
DO_R2[j] = 1 - sum((DO_mod_vec - DO_obs_vec) .* (DO_mod_vec - DO_obs_vec)) / sum((DO_obs_vec - mean(DO_obs_vec)) .* (DO_obs_vec - mean(DO_obs_vec)));
}
}
****