Computer experiments
We start by generating a maximin Latin hypercube, following the run of thumb of having 10 times as many runs as input variables. The function maximinSA_LHS requires a starting design for the optimisation; we will use a random Latin hypercube of the correct size. We can visualise our final design by looking at one-dimensional projections (which should be uniform) and two-dimensional projections (the design criterion doesn’t consider these, but they hopefully will look OK).
library(DiceDesign)
start.design <- lhsDesign(80, 8)$design
maximin.d <- maximinSA_LHS(start.design)
colnames(maximin.d$design) <- c("rw", "r", "Tu", "Hu", "Tl", "Hl", "L", "Kw")
hist(maximin.d$design[, 1])
pairs(maximin.d$design)
To run the experiment, we need to scale our design to the appropriate ranges (DiceDesign finds designs on the scale \([0,1]\)) and pass each row of the design to the borehole function.
library(apts.doe)
l <- c(.05, 100, 63070, 990, 63.1, 700, 1120, 9855)
u <- c(.15, 50000, 115600, 1110, 116, 820, 1680, 12045)
center <- -l / (u - l)
scalev <- 1 / (u - l)
bhdesign <- scale(maximin.d$design, center = center, scale = scalev)
bhdesign
## rw r Tu Hu Tl Hl L
## [1,] 0.11537239 29219.0282 75382.12 1068.3173 100.50402 709.3176 1215.375
## [2,] 0.14546237 33514.5571 85695.92 993.0483 104.21359 801.9650 1439.911
## [3,] 0.05535548 10359.4135 65444.77 1049.6694 71.83889 813.4733 1474.682
## [4,] 0.12662909 46364.5476 67959.53 1018.0191 83.26376 744.0609 1226.499
## [5,] 0.14826760 17918.0545 102395.51 1043.3945 88.60545 755.9380 1466.979
## [6,] 0.07673155 37013.2550 103611.68 1078.5740 94.62532 703.0235 1149.465
## [7,] 0.13140425 5431.1244 71252.47 1060.4337 104.08928 748.3192 1377.999
## [8,] 0.09267319 45991.7186 92676.63 1048.3820 86.82722 753.9493 1556.455
## [9,] 0.06828096 7684.0230 82597.72 1001.8662 77.86824 793.1268 1199.162
## [10,] 0.10024538 17454.0793 106630.97 1029.5972 74.62280 791.0753 1610.302
## [11,] 0.10745694 23907.8701 74227.83 1093.4344 64.56479 784.1258 1645.538
## [12,] 0.14046010 25035.2237 94832.93 1000.0361 73.79977 798.8818 1320.860
## [13,] 0.12378778 6696.8456 105347.17 1070.7515 112.68672 720.7161 1536.508
## [14,] 0.08256554 35750.2627 90150.01 1046.2768 109.84955 806.0497 1335.705
## [15,] 0.10843667 43652.1238 80056.60 1041.6872 63.89793 819.7385 1409.281
## [16,] 0.11915562 216.4491 90991.28 1007.6579 85.54829 736.6169 1444.136
## [17,] 0.08726074 32474.8338 96346.25 1087.1908 83.61669 789.0849 1171.394
## [18,] 0.05024158 41304.5375 63979.75 1060.7550 110.53604 817.4315 1597.155
## [19,] 0.10336697 4582.1316 109748.44 1051.8512 87.94302 723.3986 1551.503
## [20,] 0.11625584 13937.5755 69929.21 1081.5701 101.84764 802.1231 1155.981
## [21,] 0.09530599 38821.8122 113310.53 1015.1067 67.68323 773.3843 1249.737
## [22,] 0.11026803 48228.7805 89150.53 1078.3435 106.01184 727.9497 1384.578
## [23,] 0.12584459 34047.9232 91695.07 1098.9901 65.62162 715.6370 1271.196
## [24,] 0.05690062 12267.1576 76991.21 1004.0658 84.85474 732.5157 1414.009
## [25,] 0.05323357 18909.4550 104705.24 1025.9688 80.74662 816.8315 1431.135
## [26,] 0.12325839 47955.7351 108796.01 1106.1315 106.93599 781.4195 1206.240
## [27,] 0.13068331 22908.0898 111607.97 994.8169 81.17870 724.4976 1490.263
## [28,] 0.06568001 19807.5055 98743.60 991.1691 70.30063 712.5901 1134.984
## [29,] 0.14462486 39755.8378 78584.60 1030.8242 96.43516 788.3104 1310.067
## [30,] 0.08135119 9676.3349 114436.25 1026.7778 85.88069 758.4766 1306.977
## [31,] 0.14883191 11695.1538 85029.02 1085.3877 79.69802 779.8023 1164.253
## [32,] 0.13807450 31158.4296 67064.50 1108.3249 113.38561 747.2191 1391.777
## [33,] 0.07036292 32882.4044 87082.03 1033.3715 63.58912 795.9810 1584.049
## [34,] 0.14740286 15756.0164 80528.32 1067.0936 90.90622 800.2501 1502.210
## [35,] 0.08556261 49572.4357 113804.02 1034.4461 93.15134 778.3113 1520.349
## [36,] 0.09918574 26619.5901 63528.95 998.1053 68.73569 711.3410 1327.345
## [37,] 0.11429947 36424.4928 95442.48 1083.0503 87.48217 777.3974 1669.518
## [38,] 0.10468956 45326.4390 79146.04 1062.2473 69.07220 707.6860 1513.953
## [39,] 0.05499937 28357.9052 101749.60 1091.2152 89.46995 716.5142 1652.387
## [40,] 0.07286393 31450.7610 74278.94 1075.9705 79.04955 719.2500 1276.114
## [41,] 0.09825578 20667.2211 89889.52 1039.4895 65.88083 737.7054 1219.017
## [42,] 0.09208811 40974.9448 69006.64 1024.3730 106.27045 729.3686 1511.314
## [43,] 0.09656920 12720.2762 100994.63 1005.1530 77.32677 811.6852 1405.582
## [44,] 0.06262735 27124.6050 99465.08 1037.0555 110.73302 746.4228 1143.578
## [45,] 0.08110154 8460.7936 72227.38 1095.0119 114.06277 721.2008 1345.387
## [46,] 0.14297965 35396.1110 77732.18 1053.4431 70.57050 765.2836 1295.570
## [47,] 0.08953049 40625.6981 83868.31 1009.8812 115.01395 700.5589 1179.720
## [48,] 0.11275790 1422.6345 88657.36 1028.2201 73.47946 754.8292 1617.818
## [49,] 0.12888255 41903.5233 87806.97 1019.3981 75.74385 706.8625 1125.641
## [50,] 0.09035706 44222.8789 64845.37 1081.4846 89.78920 770.4595 1364.803
## [51,] 0.05919128 15646.4149 102482.22 1103.1121 92.70216 713.8839 1262.364
## [52,] 0.06175129 16678.5523 94089.25 1016.3615 90.49860 733.3580 1572.307
## [53,] 0.07542878 13607.2687 83226.52 1088.1101 94.15689 759.0658 1609.237
## [54,] 0.10884000 25414.8398 66545.16 1099.7812 103.06002 702.5487 1592.742
## [55,] 0.12197877 38658.5840 104181.05 1035.3087 102.51564 743.3940 1342.045
## [56,] 0.13619583 21272.8531 115399.87 1097.5579 75.64264 796.6270 1254.332
## [57,] 0.05185551 27572.1276 81283.59 1105.1241 105.29440 749.6732 1626.149
## [58,] 0.11786273 23287.1938 105979.37 1109.8261 113.34630 791.8956 1526.425
## [59,] 0.13357524 18555.2707 107552.39 1093.7044 71.14706 751.2677 1395.319
## [60,] 0.06107568 11012.2957 86576.20 1008.5761 111.43871 806.5267 1282.055
## [61,] 0.05844912 42783.6494 98336.54 1057.8017 78.33010 730.8198 1366.445
## [62,] 0.12122857 21859.8515 109093.34 1012.3711 108.34387 771.4973 1635.609
## [63,] 0.07837248 14477.9223 70794.39 1056.3205 82.05365 768.8164 1566.171
## [64,] 0.09381159 34880.8278 107959.41 1051.4048 98.96874 775.3314 1288.599
## [65,] 0.13904817 9375.8857 93330.19 1054.5696 95.06026 739.8963 1239.122
## [66,] 0.07926744 26056.8854 81759.95 1040.9956 98.75985 815.1655 1453.896
## [67,] 0.10155085 2274.8016 84508.31 1066.1041 108.93848 766.4887 1661.411
## [68,] 0.10530698 1150.7812 97625.51 1071.7642 98.10596 803.6990 1427.224
## [69,] 0.06930397 7492.2642 97099.05 1090.2116 101.35303 785.7663 1185.073
## [70,] 0.14191578 49182.9315 100272.85 1064.5311 96.88066 810.4880 1478.020
## [71,] 0.08811211 6277.1382 76088.29 1101.0572 66.98848 774.3849 1190.959
## [72,] 0.06674825 29765.0298 68383.32 1044.1148 115.68107 741.8344 1545.613
## [73,] 0.13473304 22164.6973 65854.30 1073.8928 68.32198 782.6980 1233.671
## [74,] 0.11179465 38069.6631 76504.51 996.9301 72.45594 760.2474 1494.353
## [75,] 0.07217883 30415.7865 112836.86 992.6525 82.30438 763.2768 1641.383
## [76,] 0.08425932 3941.1619 73513.91 1020.9802 95.76458 705.5891 1129.240
## [77,] 0.07486668 3695.3500 112285.99 1074.1623 76.39219 761.5849 1357.239
## [78,] 0.06421231 46931.9994 92319.58 1013.7209 107.91794 735.9758 1575.166
## [79,] 0.12818286 44915.9181 72855.72 1021.5587 99.57134 726.1911 1678.686
## [80,] 0.13744039 3093.0436 110982.75 1002.3575 92.08545 808.8613 1460.183
## Kw
## [1,] 12018.159
## [2,] 10659.358
## [3,] 10371.356
## [4,] 11097.064
## [5,] 10159.971
## [6,] 10768.605
## [7,] 10238.328
## [8,] 9975.980
## [9,] 11793.153
## [10,] 11848.767
## [11,] 10377.094
## [12,] 10969.724
## [13,] 11927.574
## [14,] 11592.299
## [15,] 11478.995
## [16,] 12007.658
## [17,] 11886.046
## [18,] 10451.922
## [19,] 10874.200
## [20,] 11440.417
## [21,] 11055.051
## [22,] 10721.943
## [23,] 11712.842
## [24,] 10423.417
## [25,] 11028.423
## [26,] 10612.965
## [27,] 11189.162
## [28,] 11547.664
## [29,] 11936.648
## [30,] 11407.783
## [31,] 10472.416
## [32,] 11314.454
## [33,] 10581.374
## [34,] 11062.217
## [35,] 11164.030
## [36,] 10993.274
## [37,] 11990.246
## [38,] 11234.917
## [39,] 10691.952
## [40,] 10303.823
## [41,] 10548.950
## [42,] 11869.202
## [43,] 9867.160
## [44,] 11466.005
## [45,] 11751.526
## [46,] 10087.467
## [47,] 10931.112
## [48,] 10043.607
## [49,] 10140.466
## [50,] 11273.097
## [51,] 11673.577
## [52,] 11739.460
## [53,] 10115.595
## [54,] 10634.730
## [55,] 11634.032
## [56,] 9884.088
## [57,] 11572.273
## [58,] 10270.492
## [59,] 10900.585
## [60,] 10866.570
## [61,] 11809.990
## [62,] 11215.913
## [63,] 11648.332
## [64,] 10007.112
## [65,] 11377.770
## [66,] 9945.304
## [67,] 11347.897
## [68,] 9919.965
## [69,] 10798.987
## [70,] 11283.628
## [71,] 10061.167
## [72,] 10198.403
## [73,] 11509.160
## [74,] 10332.810
## [75,] 10487.524
## [76,] 10817.402
## [77,] 10229.048
## [78,] 10745.249
## [79,] 10521.850
## [80,] 11128.587
## attr(,"scaled:center")
## [1] -0.500000000 -0.002004008 -1.200647249 -8.250000000 -1.192816635
## [6] -5.833333333 -2.000000000 -4.500000000
## attr(,"scaled:scale")
## [1] 1.000000e+01 2.004008e-05 1.903674e-05 8.333333e-03 1.890359e-02
## [6] 8.333333e-03 1.785714e-03 4.566210e-04
y <- apply(bhdesign, 1, borehole)
y
## [1] 147.24706 93.16504 15.96210 123.48550 136.12401 64.82985 124.97501
## [8] 50.70603 29.96537 55.07848 70.27019 102.17133 130.05917 44.49503
## [15] 66.14203 100.00442 71.93234 12.61441 76.92288 116.53137 60.55221
## [22] 103.05302 173.32664 20.32551 14.32303 135.26979 107.86412 38.26139
## [29] 143.42115 48.49855 188.62181 174.35139 24.58111 132.77665 43.09265
## [36] 72.86306 89.48699 89.90764 23.00756 47.87238 78.59548 61.47283
## [43] 39.59209 35.82351 67.27342 142.10412 71.87257 67.45110 145.10412
## [50] 65.57051 39.53133 25.26330 36.88517 98.24134 117.35111 136.60791
## [57] 21.34307 92.89031 148.23353 20.03475 30.25394 75.83410 41.12257
## [64] 59.01433 173.69731 30.40927 66.08662 64.67514 41.75676 121.48120
## [71] 66.94430 27.87099 152.65194 63.77259 23.92452 67.12967 41.35870
## [78] 24.49958 94.93653 86.83445
We use DiceKriging and the km function to fit a GP model to this data. We use our original design on the interval \([0,1]\) so that the correlation parameters are comparable (they are all on the same scale).
library(DiceKriging)
gp <- km(formula = ~., design = maximin.d$design, response = y, control = list(trace = F))
gp
##
## Call:
## km(formula = ~., design = maximin.d$design, response = y, control = list(trace = F))
##
## Trend coeff.:
## Estimate
## (Intercept) 21.8785
## rw 144.1632
## r 0.8479
## Tu -1.8313
## Hu 35.0752
## Tl 1.6642
## Hl -30.1165
## L -32.6835
## Kw 14.6344
##
## Covar. type : matern5_2
## Covar. coeff.:
## Estimate
## theta(rw) 0.4887
## theta(r) 1.9782
## theta(Tu) 1.9749
## theta(Hu) 1.4027
## theta(Tl) 1.9694
## theta(Hl) 1.4490
## theta(L) 1.1735
## theta(Kw) 1.9644
##
## Variance estimate: 163.1625
The dominant variable is \(r_w\), with a large linear coefficient and a small \(\theta\) (correlation parameter).
We can graphically explore a bit using leave-one-out predictions. If we plot these against the response, we would hope to see a line close to the diagonal, with narrow intervals mostly crossing the line (so the response is usually within the interval). We can also work out an approximate \(R^2\) value as the square of the correlation between the response and the mean prediction. The model looks pretty good on these measures.
library(ggplot2)
loo <- leaveOneOut.km(gp, type = "UK", trend.reestim = T)
summary_df <- data.frame(y = y, pred = loo$mean, ymin = loo$mean - 2 * loo$sd, ymax = loo$mean + 2 * loo$sd)
ggplot(summary_df, aes(x = y, y = pred, color = 100 * maximin.d$design[, "rw"])) + geom_point(size = 2) +
geom_abline(slope = 1, intercept = 0) +
geom_errorbar(aes(ymin = ymin, ymax = ymax)) + labs(color = "r_w")
with(summary_df, cor(y, pred)) ^ 2
## [1] 0.9979813