buildsim calibration

Integrate BC-STAN scripts from AChong.

Author: weilixu

.idea/workspace.xml

@@ -0,0 +1,121 @@
+yc,xc1,xc2,xc3,tc1,tc2,tc3
+1836.311382,1.691666667,37.93010753,204.0215054,0.726,19.535,0.0011998
+1600.207111,0.061011905,51.19196429,181.8467262,0.726,19.535,0.0011998
+1738.757931,4.350672043,50.60752688,219.4072581,0.726,19.535,0.0011998
+1682.950916,5.442638889,56.03472222,199.3166667,0.726,19.535,0.0011998
+1775.641414,14.06424731,44.78091398,244.1129032,0.726,19.535,0.0011998
+2024.757846,21.35458333,37.57916667,305.2166667,0.726,19.535,0.0011998
+2105.986898,22.26639785,49.41532258,272.4677419,0.726,19.535,0.0011998
+2138.781979,22.61505376,42.75537634,222.2553763,0.726,19.535,0.0011998
+1832.845423,18.45805556,53.5,225.0458333,0.726,19.535,0.0011998
+1652.636921,8.881451613,54.77016129,233.7217742,0.726,19.535,0.0011998
+1641.316706,6.096666667,42.07638889,217.4472222,0.726,19.535,0.0011998
+1915.45202,0.197043011,48.17607527,224.5537634,0.726,19.535,0.0011998
+2106.973629,1.691666667,37.93010753,204.0215054,0.58,29.558,0.0006871
+1824.393576,0.061011905,51.19196429,181.8467262,0.58,29.558,0.0006871
+2038.370765,4.350672043,50.60752688,219.4072581,0.58,29.558,0.0006871
+1952.372714,5.442638889,56.03472222,199.3166667,0.58,29.558,0.0006871
+2130.126984,14.06424731,44.78091398,244.1129032,0.58,29.558,0.0006871
+2397.227572,21.35458333,37.57916667,305.2166667,0.58,29.558,0.0006871
+2489.730744,22.26639785,49.41532258,272.4677419,0.58,29.558,0.0006871
+2533.991209,22.61505376,42.75537634,222.2553763,0.58,29.558,0.0006871
+2189.037912,18.45805556,53.5,225.0458333,0.58,29.558,0.0006871
+1989.524663,8.881451613,54.77016129,233.7217742,0.58,29.558,0.0006871
+1943.786733,6.096666667,42.07638889,217.4472222,0.58,29.558,0.0006871
+2182.146452,0.197043011,48.17607527,224.5537634,0.58,29.558,0.0006871
+1976.665565,1.691666667,37.93010753,204.0215054,0.422,24.97,0.0001627
+1706.465065,0.061011905,51.19196429,181.8467262,0.422,24.97,0.0001627
+1926.973194,4.350672043,50.60752688,219.4072581,0.422,24.97,0.0001627
+1839.710073,5.442638889,56.03472222,199.3166667,0.422,24.97,0.0001627
+2028.951466,14.06424731,44.78091398,244.1129032,0.422,24.97,0.0001627
+2271.219777,21.35458333,37.57916667,305.2166667,0.422,24.97,0.0001627
+2361.691578,22.26639785,49.41532258,272.4677419,0.422,24.97,0.0001627
+2398.918526,22.61505376,42.75537634,222.2553763,0.422,24.97,0.0001627
+2081.633159,18.45805556,53.5,225.0458333,0.422,24.97,0.0001627
+1899.969234,8.881451613,54.77016129,233.7217742,0.422,24.97,0.0001627
+1847.183727,6.096666667,42.07638889,217.4472222,0.422,24.97,0.0001627
+2036.259799,0.197043011,48.17607527,224.5537634,0.422,24.97,0.0001627
+1645.848066,1.691666667,37.93010753,204.0215054,0.706,12.971,0.0011428
+1432.533047,0.061011905,51.19196429,181.8467262,0.706,12.971,0.0011428
+1544.25852,4.350672043,50.60752688,219.4072581,0.706,12.971,0.0011428
+1499.753292,5.442638889,56.03472222,199.3166667,0.706,12.971,0.0011428
+1561.758979,14.06424731,44.78091398,244.1129032,0.706,12.971,0.0011428
+1791.965459,21.35458333,37.57916667,305.2166667,0.706,12.971,0.0011428
+1866.808322,22.26639785,49.41532258,272.4677419,0.706,12.971,0.0011428
+1891.104205,22.61505376,42.75537634,222.2553763,0.706,12.971,0.0011428
+1617.316513,18.45805556,53.5,225.0458333,0.706,12.971,0.0011428
+1453.964856,8.881451613,54.77016129,233.7217742,0.706,12.971,0.0011428
+1455.169041,6.096666667,42.07638889,217.4472222,0.706,12.971,0.0011428
+1732.329969,0.197043011,48.17607527,224.5537634,0.706,12.971,0.0011428
+1473.378259,1.691666667,37.93010753,204.0215054,0.413,6.012,0.0002739
+1248.091943,0.061011905,51.19196429,181.8467262,0.413,6.012,0.0002739
+1381.207432,4.350672043,50.60752688,219.4072581,0.413,6.012,0.0002739
+1335.370215,5.442638889,56.03472222,199.3166667,0.413,6.012,0.0002739
+1425.239521,14.06424731,44.78091398,244.1129032,0.413,6.012,0.0002739
+1626.465221,21.35458333,37.57916667,305.2166667,0.413,6.012,0.0002739
+1697.916461,22.26639785,49.41532258,272.4677419,0.413,6.012,0.0002739
+1712.343273,22.61505376,42.75537634,222.2553763,0.413,6.012,0.0002739
+1478.975722,18.45805556,53.5,225.0458333,0.413,6.012,0.0002739
+1343.705986,8.881451613,54.77016129,233.7217742,0.413,6.012,0.0002739
+1318.136957,6.096666667,42.07638889,217.4472222,0.413,6.012,0.0002739
+1561.146862,0.197043011,48.17607527,224.5537634,0.413,6.012,0.0002739
+2975.979622,1.691666667,37.93010753,204.0215054,0.128,8.231,0.000735
+2600.833327,0.061011905,51.19196429,181.8467262,0.128,8.231,0.000735
+2927.14945,4.350672043,50.60752688,219.4072581,0.128,8.231,0.000735
+2845.024497,5.442638889,56.03472222,199.3166667,0.128,8.231,0.000735
+3138.696984,14.06424731,44.78091398,244.1129032,0.128,8.231,0.000735
+3504.795597,21.35458333,37.57916667,305.2166667,0.128,8.231,0.000735
+3686.821406,22.26639785,49.41532258,272.4677419,0.128,8.231,0.000735
+3707.505279,22.61505376,42.75537634,222.2553763,0.128,8.231,0.000735
+3237.556406,18.45805556,53.5,225.0458333,0.128,8.231,0.000735
+2941.936302,8.881451613,54.77016129,233.7217742,0.128,8.231,0.000735
+2818.988488,6.096666667,42.07638889,217.4472222,0.128,8.231,0.000735
+3059.892679,0.197043011,48.17607527,224.5537634,0.128,8.231,0.000735
+1835.621781,1.691666667,37.93010753,204.0215054,0.816,21.478,0.0011063
+1596.089255,0.061011905,51.19196429,181.8467262,0.816,21.478,0.0011063
+1738.643376,4.350672043,50.60752688,219.4072581,0.816,21.478,0.0011063
+1680.349216,5.442638889,56.03472222,199.3166667,0.816,21.478,0.0011063
+1779.54957,14.06424731,44.78091398,244.1129032,0.816,21.478,0.0011063
+2025.88739,21.35458333,37.57916667,305.2166667,0.816,21.478,0.0011063
+2106.001866,22.26639785,49.41532258,272.4677419,0.816,21.478,0.0011063
+2140.918137,22.61505376,42.75537634,222.2553763,0.816,21.478,0.0011063
+1835.02372,18.45805556,53.5,225.0458333,0.816,21.478,0.0011063
+1655.86323,8.881451613,54.77016129,233.7217742,0.816,21.478,0.0011063
+1642.213572,6.096666667,42.07638889,217.4472222,0.816,21.478,0.0011063
+1911.073946,0.197043011,48.17607527,224.5537634,0.816,21.478,0.0011063
+2579.296213,1.691666667,37.93010753,204.0215054,0.157,13.803,0.0004071
+2247.654682,0.061011905,51.19196429,181.8467262,0.157,13.803,0.0004071
+2533.791274,4.350672043,50.60752688,219.4072581,0.157,13.803,0.0004071
+2451.233254,5.442638889,56.03472222,199.3166667,0.157,13.803,0.0004071
+2706.172611,14.06424731,44.78091398,244.1129032,0.157,13.803,0.0004071
+3019.858217,21.35458333,37.57916667,305.2166667,0.157,13.803,0.0004071
+3166.837609,22.26639785,49.41532258,272.4677419,0.157,13.803,0.0004071
+3192.20198,22.61505376,42.75537634,222.2553763,0.157,13.803,0.0004071
+2785.337744,18.45805556,53.5,225.0458333,0.157,13.803,0.0004071
+2537.351452,8.881451613,54.77016129,233.7217742,0.157,13.803,0.0004071
+2439.359501,6.096666667,42.07638889,217.4472222,0.157,13.803,0.0004071
+2649.309258,0.197043011,48.17607527,224.5537634,0.157,13.803,0.0004071
+1584.086649,1.691666667,37.93010753,204.0215054,0.679,9.01,0.0013563
+1390.610008,0.061011905,51.19196429,181.8467262,0.679,9.01,0.0013563
+1484.734307,4.350672043,50.60752688,219.4072581,0.679,9.01,0.0013563
+1442.803491,5.442638889,56.03472222,199.3166667,0.679,9.01,0.0013563
+1481.777434,14.06424731,44.78091398,244.1129032,0.679,9.01,0.0013563
+1712.476965,21.35458333,37.57916667,305.2166667,0.679,9.01,0.0013563
+1786.29391,22.26639785,49.41532258,272.4677419,0.679,9.01,0.0013563
+1805.557701,22.61505376,42.75537634,222.2553763,0.679,9.01,0.0013563
+1539.256742,18.45805556,53.5,225.0458333,0.679,9.01,0.0013563
+1378.373847,8.881451613,54.77016129,233.7217742,0.679,9.01,0.0013563
+1386.676202,6.096666667,42.07638889,217.4472222,0.679,9.01,0.0013563
+1680.499585,0.197043011,48.17607527,224.5537634,0.679,9.01,0.0013563
+2173.918568,1.691666667,37.93010753,204.0215054,0.384,23.152,0.0008346
+1881.648491,0.061011905,51.19196429,181.8467262,0.384,23.152,0.0008346
+2097.010598,4.350672043,50.60752688,219.4072581,0.384,23.152,0.0008346
+2016.711717,5.442638889,56.03472222,199.3166667,0.384,23.152,0.0008346
+2191.979488,14.06424731,44.78091398,244.1129032,0.384,23.152,0.0008346
+2472.120012,21.35458333,37.57916667,305.2166667,0.384,23.152,0.0008346
+2573.891277,22.26639785,49.41532258,272.4677419,0.384,23.152,0.0008346
+2611.76857,22.61505376,42.75537634,222.2553763,0.384,23.152,0.0008346
+2258.167906,18.45805556,53.5,225.0458333,0.384,23.152,0.0008346
+2051.308275,8.881451613,54.77016129,233.7217742,0.384,23.152,0.0008346
+2001.234627,6.096666667,42.07638889,217.4472222,0.384,23.152,0.0008346
+2252.17893,0.197043011,48.17607527,224.5537634,0.384,23.152,0.0008346

buildsimcalibration/DATACOMP.csv

@@ -0,0 +1,13 @@
+yf,xf1,xf2,xf3
+1562,1.69166666666667,37.9301075268817,204.021505376344
+1395,0.0610119047619034,51.1919642857143,181.846726190476
+1279,4.35067204301076,50.6075268817204,219.407258064516
+1347,5.44263888888889,56.0347222222222,199.316666666667
+1363,14.064247311828,44.7809139784946,244.112903225806
+1691,21.3545833333333,37.5791666666667,305.216666666667
+1852,22.2663978494624,49.4153225806452,272.467741935484
+1983,22.6150537634409,42.755376344086,222.255376344086
+1650,18.4580555555556,53.5,225.045833333333
+1385,8.88145161290323,54.7701612903226,233.721774193548
+1347,6.09666666666667,42.0763888888889,217.447222222222
+1524,0.197043010752687,48.1760752688172,224.55376344086

buildsimcalibration/DATAFIELD.csv

@@ -0,0 +1,8 @@
+"""
+This package is developed jointly between BuildSim team and Professor Adrian Chong.
+
+The original script and example can be found at GITHUB page: https://github.com/adChong/bc-stan
+
+The paper published can be found at: https://www.sciencedirect.com/science/article/pii/S0378778818307539?via%3Dihub
+
+"""

buildsimcalibration/__init__.py

@@ -0,0 +1,131 @@
+data {
+  int<lower=1> n; // number of field data
+  int<lower=1> m; // number of computer simulation
+  int<lower=1> n_pred; // number of predictions
+  int<lower=1> p; // number of observable inputs x
+  int<lower=1> q; // number of calibration parameters t
+  vector[n] y; // field observations
+  vector[m] eta; // output of computer simulations
+  matrix[n, p] xf; // observable inputs corresponding to y
+  // (xc, tc): design points corresponding to eta
+  matrix[m, p] xc; 
+  matrix[m, q] tc; 
+  // x_pred: new design points for predictions
+  matrix[n_pred, p] x_pred; 
+}
+
+transformed data {
+  int<lower = 1> N;
+  vector[n+m] y_eta;
+  vector[n+m+n_pred] mu; // mean vector
+  matrix[n+n_pred, p] X; // X=[xf, x_pred]
+  
+  N = n + m + n_pred;
+  // set mean vector to zero
+  for (i in 1:N) {
+    mu[i] = 0;
+  }
+  X = append_row(xf, x_pred);
+  y_eta = append_row(y, eta); // y_eta = [y, eta]
+}
+
+parameters {
+  // tf: calibration parameters
+  // rho_eta: reparameterization of beta_eta
+  // rho_delta: reparameterization of beta_delta
+  // lambda_eta: precision parameter for eta
+  // lambda_delta: precision parameter for bias term
+  // lambda_e: precision parameter of observation error
+  // y_pred: predictions
+  row_vector<lower=0, upper=1>[q] tf; 
+  row_vector<lower=0, upper=1>[p+q] rho_eta; 
+  row_vector<lower=0, upper=1>[p] rho_delta; 
+  real<lower=0> lambda_eta; 
+  real<lower=0> lambda_delta;
+  real<lower=0> lambda_e; 
+  vector[n_pred] y_pred; 
+}
+
+transformed parameters {
+  // beta_delta: correlation parameter for bias term
+  // beta_e: correlation parameter of observation error
+  row_vector[p+q] beta_eta;
+  row_vector[p] beta_delta;
+  beta_eta = -4.0 * log(rho_eta);
+  beta_delta = -4.0 * log(rho_delta);
+}
+
+model {
+  // declare variables
+  matrix[N, p+q] xt; 
+  matrix[N, N] sigma_eta; // simulator covarinace
+  matrix[n+n_pred, n+n_pred] sigma_delta; // bias term covariance
+  matrix[N, N] sigma_z; // covariance matrix
+  matrix[N, N] L; // cholesky decomposition of covariance matrix 
+  vector[N] z; // z = [y, eta, y_pred]
+  row_vector[p] temp_delta;
+  row_vector[p+q] temp_eta;
+
+  z = append_row(y_eta, y_pred); // z = [y, eta, y_pred]
+
+  // xt = [[xf,tf],[xc,tc],[x_pred,tf]]
+  xt[1:n, 1:p] = xf;
+  xt[1:n, (p+1):(p+q)] = rep_matrix(tf, n);
+  xt[(n+1):(n+m), 1:p] = xc;
+  xt[(n+1):(n+m), (p+1):(p+q)] = tc;
+  xt[(n+m+1):N, 1:p] = x_pred;
+  xt[(n+m+1):N, (p+1):(p+q)] = rep_matrix(tf, n_pred);
+  
+  // 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_delta), 
+    n+n_pred));
+  
+  // off-diagonal elements of sigma_delta
+  for (i in 1:(n+n_pred-1)) {
+    for (j in (i+1):(n+n_pred)) {
+      temp_delta = X[i] - X[j];
+      sigma_delta[i, j] = beta_delta .* temp_delta * temp_delta';
+      sigma_delta[i, j] = exp(-sigma_delta[i, j]) / lambda_delta;
+      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[1:n, 1:n];
+  sigma_z[1:n, (n+m+1):N] = sigma_eta[1:n, (n+m+1):N] + 
+    sigma_delta[1:n, (n+1):(n+n_pred)];
+  sigma_z[(n+m+1):N, 1:n] = sigma_eta[(n+m+1):N, 1:n] + 
+    sigma_delta[(n+1):(n+n_pred),1:n];
+  sigma_z[(n+m+1):N, (n+m+1):N] = sigma_eta[(n+m+1):N, (n+m+1):N] + 
+    sigma_delta[(n+1):(n+n_pred), (n+1):(n+n_pred)];
+
+  // add observation errors
+  for (i in 1:n) {
+    sigma_z[i, i] = sigma_z[i, i] + (1.0 / lambda_e);
+  }  
+
+  // Specify priors here
+  rho_eta[1:(p+q)] ~ beta(1.0, 0.3);
+  rho_delta[1:p] ~ beta(1.0, 0.3);
+  lambda_eta ~ gamma(10, 10); // gamma (shape, rate)
+  lambda_delta ~ gamma(10, 0.3); 
+  lambda_e ~ gamma(10, 0.03); 
+
+  L = cholesky_decompose(sigma_z); // cholesky decomposition 
+  z ~ multi_normal_cholesky(mu, L);
+}