Hi, thanks for responding. I’m not quite sure which bit of the code is most useful? Here is the Stan file -
data {
int<lower=0> n; // number of field data
int<lower=0> m; // number of computer simulation
int<lower=0> p; // number of observable inputs x
int<lower=0> q; // number of calibration parameters t
vector[n] y; // field observations
vector[m] eta; // output of computer simulations
vector[n] xf; // observable inputs corresponding to y
// (xc, tc): design points corresponding to eta
vector[m] xc;
matrix[m, q] tc;
// inputs for prior distributions !!NEW!!
row_vector[q] mu_theta; // mean values for theta priors
row_vector[q] sd_theta; // standard deviations for theta priors
real shape_eta; // shape parameter for prior of emulator
real rate_eta; // rate parameter for prior of emulator
real shape_b; // shape parameter for prior of model bias
real rate_b; // rate parameter for prior of model bias
real shape_e; // shape parameter for prior of random error
real rate_e; // rate parameter for prior of random error
row_vector[p+q] a_eta; // a values for beta_eta priors
row_vector[p+q] b_eta; // b values for beta_eta priors
real a_bias; // a value for beta_b prior
real b_bias; // b value for beta_b prior
}
transformed data {
int<lower=1> N;
vector[n+m] z; // z = [y, eta]
vector[n+m] mu; // mean vector
N = n + m;
// set mean vector to zero
for (i in 1:N) {
mu[i] = 0;
}
z = append_row(y, eta);
}
parameters {
// theta: calibration parameters
// rho_eta: reparameterization of beta_eta
// rho_b: reparameterization of beta_b
// lambda_eta: precision parameter for eta
// lambda_b: precision parameter for bias term
// lambda_e: precision parameter of observation error
row_vector<lower=0,upper=1>[p+q] rho_eta;
real<lower=0,upper=1> rho_b;
real<lower=0> lambda_eta;
real<lower=0> lambda_e;
real<lower=0> lambda_b;
row_vector<lower=0,upper=1>[q] theta;
}
transformed parameters {
// beta_delta: correlation parameter for bias term
// beta_e: correlation parameter of observation error
row_vector[p+q] beta_eta;
real beta_delta;
beta_eta = -4.0 * log2(rho_eta); // !!NEW!! log2 instead of log
beta_delta = -4.0 * log2(rho_b); // !!NEW!! log2 instead of log
}
model {
// declare variables
matrix[N, (p+q)] xt;
matrix[N, N] sigma_eta; // simulator covariance
matrix[n, n] sigma_delta; // bias term covariance
matrix[N, N] sigma_z; // covariance matrix
matrix[N, N] L; // cholesky decomposition of covariance matrix
real temp_delta;
row_vector[p+q] temp_eta;
// xt = [[xt,theta],[xc,tc]]
xt[1:n, 1] = xf;
xt[(n+1):N, 1] = xc;
xt[1:n, (p+1):(p+q)] = rep_matrix(theta, n);
xt[(n+1):N, (p+1):(p+q)] = tc;
// diagonal elements of sigma_eta
sigma_eta = diag_matrix(rep_vector((1 / lambda_eta), N));
// off-diagonal elements of sigma_eta
for (i in 1:(N-1)) {
for (j in (i+1):N) {
temp_eta = xt[i] - xt[j];
sigma_eta[i, j] = beta_eta .* temp_eta * temp_eta’;
sigma_eta[i, j] = exp(-sigma_eta[i, j]) / lambda_eta;
sigma_eta[j, i] = sigma_eta[i, j];
}
}
// diagonal elements of sigma_delta
sigma_delta = diag_matrix(rep_vector((1 / lambda_b), n));
// off-diagonal elements of sigma_delta
for (i in 1:(n-1)) {
for (j in (i+1):n) {
temp_delta = xf[i] - xf[j];
sigma_delta[i, j] = beta_delta* temp_delta * temp_delta’;
sigma_delta[i, j] = exp(-sigma_delta[i, j]) / lambda_b;
sigma_delta[j, i] = sigma_delta[i, j];
}
}
// computation of covariance matrix sigma_z
sigma_z = sigma_eta;
sigma_z[1:n, 1:n] = sigma_eta[1:n, 1:n] + sigma_delta;
// add observation errors
for (i in 1:n) {
sigma_z[i, i] = sigma_z[i, i] + (1.0 / lambda_e);
}
// Specify distributions for thetas and hyperparameters !!NEW!!
theta[1:q] ~ normal(mu_theta, sd_theta);
rho_eta[1:(p+q)] ~ beta(a_eta, b_eta);
rho_b ~ beta(a_bias, b_bias);
lambda_b ~ gamma(shape_b, rate_b);
lambda_eta ~ gamma(shape_eta, rate_eta);
lambda_e ~ gamma(shape_e, rate_e);
L = cholesky_decompose(sigma_z); // cholesky decomposition
z ~ multi_normal_cholesky(mu, L);
}
This is called in MATLAB by:
data = struct (‘xf’,xxf,‘xc’, xxc,‘x’,xx,‘y’,yyf,‘eta’,yyc,‘tc’,tt,… % normalized data
‘n’,n,‘m’,m,‘p’,p,‘q’,q,… % dimensions of data
‘mu_theta’,mu_theta,‘sd_theta’,sd_theta,… % parameters for priors calibration parameters
‘shape_eta’, shape_eta,‘rate_eta’,rate_eta,… % parameters for pripor hyperparameters emulator
‘a_eta’,a_eta,‘b_eta’,b_eta,…
‘shape_b’,shape_b,‘rate_b’,rate_b,… % parameters for pripor hyperparameters model bias
‘a_bias’,a_bias,‘b_bias’,b_bias,…
‘shape_e’,shape_e,‘rate_e’,rate_e); % parameters for pripor hyperparameters random error
calibration = stan (‘file’,name_stan,‘data’,data,‘iter’,iter,‘chains’,chains,‘verbose’,true,‘refresh’,(iter/100),…
‘inc_warmup’,true);
:
Any advice gratefully received, thankyou :)