The first part of the context is not the important part, but here it is:
So I'm trying to fit a matern model with the variofit function in R for a school project. When trying fix.kappa == FALSE for it to find itself the best value for kappa failed for some reason. Therefore I tried to find it manually, by trying a series of values for kappa, and searching manually for one that would minimise the MSE.
Here is the strange part : I made this code that calculates the the MSE for each value of kappa between 0.5 and 50:
kappa_list = seq(0.5, 50, length.out = 50)
MSE_list = c()
for(kappa_temp in kappa_list){
variofit_temp = variofit(
emp_variog2,
fix.nugget = FALSE, nugget = nugget_0,
fix.kappa = TRUE, kappa = kappa_temp,
ini.cov.pars = c(part_sill_0,range_0),
cov.model = "matern",
messages = F
)
MSE_list = append(variofit_temp$value, MSE_list)
}
plot(kappa_list, MSE_list)
Then I remembered that kappa could take a tinier value than 0.5, so I re wrote the code:
kappa_list = seq(0.1, 50, length.out = 50)
MSE_list = c()
for(kappa_temp in kappa_list){
variofit_temp = variofit(
emp_variog2,
fix.nugget = FALSE, nugget = nugget_0,
fix.kappa = TRUE, kappa = kappa_temp,
ini.cov.pars = c(part_sill_0,range_0),
cov.model = "matern",
messages = F
)
MSE_list = append(variofit_temp$value, MSE_list)
}
plot(kappa_list, MSE_list)
But I realised the value it shows for kappa = 50 is different for this loop than the previous one. I thought at first that it was just a precision error, because the MSE is very sensitive around 50, so I looked if the last value of the sequence was really exactly 50 or 50.00000001, but it seems that it's exactly 50.
I tried some other things, but figured out that "simulating" the loop by playing each line one by one, strangly didn't the problem ! That's the weird part. Here's the code that a played in the console :
> min_kappa = 0.1
> kappa_list = seq(min_kappa, 50, length.out = 5)
> kappa_temp = kappa_list[length(kappa_list)]
> kappa_temp
[1] 50
> variofit_temp = variofit(
+ emp_variog2,
+ fix.nugget = FALSE, nugget = nugget_0,
+ fix.kappa = TRUE, kappa = kappa_temp,
+ ini.cov.pars = c(part_sill_0,range_0),
+ cov.model = "matern",
+ messages = F
+ )
> variofit_temp$value
[1] 1721348178
> variofit_temp = variofit(
+ emp_variog2,
+ fix.nugget = FALSE, nugget = nugget_0,
+ fix.kappa = TRUE, kappa = 50,
+ ini.cov.pars = c(part_sill_0,range_0),
+ cov.model = "matern",
+ messages = F
+ )
> variofit_temp$value
[1] 1721348178
> kappa_list = seq(min_kappa, 50, length.out = 5)
> print(paste0("kappa_list = ", kappa_list))
[1] "kappa_list = 0.1" "kappa_list = 12.575" "kappa_list = 25.05" "kappa_list = 37.525" "kappa_list = 50"
> MSE_list = c()
> for(kappa_temp in kappa_list){
+ set.seed(123)
+ print(paste0("kappa_temp = ", format(kappa_temp, digits = 22, scientific = FALSE)))
+
+ kappa_temp <- round(kappa_temp, digits = 1)
+
+ variofit_temp = variofit(
+ emp_variog2,
+ fix.nugget = FALSE, nugget = nugget_0,
+ fix.kappa = TRUE, kappa = kappa_temp,
+ ini.cov.pars = c(part_sill_0,range_0),
+ cov.model = "matern",
+ messages = F
+ )
+
+ MSE_list = append(variofit_temp$value, MSE_list)
+ }
[1] "kappa_temp = 0.1000000000000000055511"
[1] "kappa_temp = 12.57499999999999928946"
[1] "kappa_temp = 25.05000000000000071054"
[1] "kappa_temp = 37.52499999999999857891"
[1] "kappa_temp = 50"
>
> print(paste0("MSE for kappa = ", kappa_list[length(kappa_list)],
+ " is ", MSE_list[length(MSE_list)]))
[1] "MSE for kappa = 50 is 21781768624.2714"
Here is the full code :
# ---- Libraries ----
library(ggplot2)
library(geoR)
library(sp)
# ---- Data importation ----
parana = structure(
list(
east = c(
402.95294, 501.70493, 556.32618, 573.40431,
702.42284, 668.54416, 435.8477, 434.01245, 432.46225, 162.14585,
169.47663, 717.84675, 713.34332, 726.66302, 618.46885, 661.88082,
506.71667, 561.7874, 585.63949, 531.71056, 308.97183, 306.10469,
150.12196, 155.3399, 620.00181, 475.05766, 366.87187, 479.49527,
615.72353, 599.07288, 566.56758, 484.6613, 445.54377, 498.29528,
477.94465, 353.18336, 338.52463, 267.23204, 676.87847, 538.88998,
340.96468, 203.93695, 254.26247, 719.46243, 768.50871, 721.65469,
687.75359, 599.06043, 541.89327, 449.6865, 357.51524, 222.33953,
297.32299, 184.17268, 401.98944, 365.36206, 331.68065, 294.95009,
674.31593, 518.6908, 562.79775, 399.91232, 309.69511, 308.21592,
397.59888, 206.29791, 289.38003, 246.80304, 279.55302, 231.82964,
245.59459, 410.80902, 400.01532, 405.11185, 422.47533, 457.73378,
388.47179, 299.07823, 320.4695, 317.7274, 328.1532, 341.96468,
346.05942, 371.75809, 377.09969, 395.46678, 244.79297, 232.08612,
212.36295, 250.92841, 249.11474, 262.98869, 264.83559, 277.18097,
267.48289, 287.42489, 291.21859, 296.594, 295.11083, 211.66139,
189.31521, 165.97002, 164.28177, 184.98705, 711.89977, 749.59433,
702.73299, 650.76351, 640.73071, 600.65702, 421.28296, 451.44393,
419.33262, 377.71707, 357.49565, 372.27262, 287.35834, 285.6835,
224.32153, 218.99256, 262.66451, 217.24028, 212.63015, 249.35852,
245.90862, 251.23978, 197.7439, 172.85298, 187.83331, 568.34196,
541.68923, 493.34111, 471.68756, 445.03235, 479.99208, 333.22941,
400.32814, 343.6607, 320.24316, 331.8482, 235.32685, 226.53333,
692.54456
),
north = c(
164.52841, 428.771, 445.27065, 447.04177,
272.2959, 261.6707, 286.54044, 317.90506, 288.37001, 286.29984,
334.53158, 181.56741, 214.88336, 207.27952, 127.43265, 128.81455,
212.86691, 138.89773, 212.59809, 131.62153, 181.97682, 148.69468,
162.12624, 154.86652, 112.64927, 83.6448399999998, 92.2122400000002,
461.96806, 400.63795, 456.12037, 445.22857, 423.22788, 400.9933,
426.92592, 356.79285, 446.50166, 381.74695, 373.39454, 254.18138,
308.77851, 271.01967, 309.41998, 240.19632, 177.84707, 199.12352,
209.20941, 187.56295, 210.66083, 185.11924, 196.16437, 188.1198,
154.59123, 176.26524, 211.03701, 70.3728200000003, 77.4256500000004,
110.26775, 113.43465, 185.89772, 378.94032, 362.2088, 354.61628,
375.85998, 359.22457, 443.1808, 359.3556, 370.04023, 374.89923,
345.88871, 450.35979, 445.05926, 242.10715, 339.85273, 338.04194,
262.48003, 306.9234, 299.16275, 285.24334, 329.84407, 279.96009,
256.09127, 271.01967, 313.5301, 278.71589, 247.38436, 267.83959,
297.28745, 343.23287, 311.44174, 238.28897, 339.84732, 321.6157,
312.41223, 281.21777, 255.20257, 273.991, 246.34864, 340.60144,
325.80878, 263.38789, 281.3855, 267.90686, 267.86614, 323.79426,
124.40728, 175.47115, 179.95796, 178.79807, 180.7536, 199.57391,
177.57037, 174.02255, 223.70296, 179.08881, 189.96577, 223.34289,
170.56964, 170.54282, 139.84991, 154.52067, 146.1385, 158.17924,
217.21353, 142.19763, 147.67415, 131.14844, 205.80201, 192.283,
198.18379, 114.86877, 120.51877, 98.4340199999995, 103.94229,
105.70264, 118.72419, 119.52053, 70.36, 84.566, 93.49421, 97.3436600000001,
92.0340700000003, 114.02601, 170.87504
),
data = c(306.09, 200.88,
167.07, 162.77, 163.57, 178.61, 301.54, 282.07, 319.12, 244.67,
233.31, 224.46, 206.12, 248.99, 237.87, 222.87, 263.1, 236.91,
247.01, 240.58, 304.28, 351.73, 277.92, 323.08, 253.32, 315.33,
379.94, 197.09, 199.91, 167, 182.88, 197.1, 257.5, 205.16, 224.07,
212.5, 242.08, 247.79, 187.27, 222.54, 313.6, 269.92, 321.69,
208.89, 238.62, 248.76, 193.48, 240.48, 265.56, 302.13, 335.41,
330.87, 329.49, 262.81, 365.88, 359.08, 344.59, 366.1, 201.84,
218.27, 200.38, 229.4, 235.07, 236.25, 228.82, 258.12, 232.17,
248.17, 240.66, 184.59, 165.46, 320.31, 232.8, 266.27, 301.1,
244.78, 248.57, 299, 269.49, 300.89, 311.81, 317.34, 286.93,
306.64, 349.09, 289.11, 280.49, 237.11, 269.68, 307.47, 253.57,
226.71, 299.69, 266.62, 331.46, 270.67, 321.15, 258.2, 288.85,
255.76, 264.77, 271.76, 313.63, 252.94, 246.75, 215.95, 299.71,
203.4, 209.47, 185.63, 304.3, 320.05, 273.01, 295.62, 334.95,
350.18, 291.35, 319.57, 367.32, 332.04, 332.18, 323.82, 334.22,
345.33, 345.24, 383.13, 348.67, 282.54, 291.07, 245.42, 291.63,
288.5, 314.65, 331.85, 305.43, 285.05, 359.51, 356.82, 372.56,
413.7, 412.44, 400.96, 199.37
)
), class = "data.frame", row.names = c(NA, -143L)
)
parana_borders = structure(
list(
east = c(
670.204875, 663.718688, 656.066688, 649.971438,
642.634563, 635.271688, 628.742125, 622.755875, 617.270813, 610.457,
604.42175, 600.262, 592.621, 585.444625, 581.92875, 574.871563,
567.99075, 565.6725, 560.5745, 553.86975, 546.261688, 542.435063,
537.837375, 534.376563, 529.409438, 528.1835, 520.268156, 513.789781,
509.833156, 505.105344, 499.977219, 491.988906, 486.094406, 478.869313,
474.045719, 470.925344, 473.2195, 476.876406, 476.248156, 468.888031,
462.077594, 457.88275, 453.517656, 448.93925, 440.925281, 433.288844,
425.385563, 417.613313, 409.689906, 402.012969, 395.090313, 388.905156,
381.620781, 373.674125, 365.810969, 358.028781, 350.718188, 343.347406,
336.159188, 329.237156, 321.329844, 313.663188, 305.80225, 298.281313,
290.861594, 284.227281, 277.077844, 271.504219, 263.543781, 256.520563,
248.8575, 241.656109, 235.520094, 231.700766, 226.547563, 225.8555,
219.657406, 216.403594, 215.340656, 216.788484, 213.318719, 212.552734,
206.664047, 201.807359, 199.6235, 192.151453, 187.422156, 188.983141,
184.399016, 178.437156, 174.949359, 170.373547, 163.837438, 157.839297,
155.134688, 150.051703, 144.086359, 138.815359, 138.297594, 137.987281,
144.653203, 149.781453, 148.382547, 153.388125, 159.686078, 155.62975,
153.212984, 154.687313, 157.737141, 163.658594, 160.675313, 164.211328,
172.130047, 165.728922, 158.935609, 157.750344, 159.271844, 158.399938,
155.783297, 157.487594, 158.247594, 162.589, 165.391719, 157.077922,
157.754734, 160.412781, 162.614203, 166.648797, 163.022797, 163.530688,
166.652172, 171.481313, 171.281031, 165.928656, 165.506375, 165.3565,
167.316813, 170.03325, 167.977563, 163.491438, 161.473719, 163.997516,
170.930609, 178.861531, 185.087641, 190.730766, 190.762938, 193.086672,
195.451828, 198.341609, 200.613078, 203.828609, 204.428203, 208.680922,
214.333047, 221.180953, 224.122406, 226.760719, 229.690422, 231.72975,
231.220078, 235.038641, 240.402063, 247.119406, 254.061875, 261.190891,
268.550156, 275.370188, 282.111844, 287.672, 293.985563, 301.211688,
307.937906, 315.517, 323.409656, 331.151406, 336.672156, 343.911,
351.417813, 356.621688, 362.441, 367.828219, 368.792531, 375.514781,
381.666813, 383.534469, 387.002438, 394.258938, 396.523313, 396.644688,
402.613531, 408.970906, 416.527875, 423.576406, 430.054938, 437.537875,
444.892563, 452.482813, 459.560219, 464.211406, 465.268031, 467.336875,
465.45875, 472.983781, 476.725063, 480.198531, 478.720594, 479.167094,
480.899469, 482.783188, 483.727938, 483.194875, 485.376531, 488.530688,
491.747375, 494.156156, 499.878813, 496.505469, 499.45075, 502.862656,
495.062781, 491.5965, 498.709813, 504.852094, 500.144813, 500.993125,
505.134438, 503.235969, 506.228938, 508.160219, 509.2235, 503.473406,
509.178344, 516.034031, 520.147219, 521.236875, 527.753063, 531.995563,
538.783625, 545.707125, 553.139188, 560.024625, 566.74675, 574.19725,
582.02375, 589.887375, 597.229313, 603.811188, 605.871, 612.286938,
613.857063, 621.287813, 628.639063, 631.817438, 637.492438, 634.20475,
629.226563, 630.906875, 638.968875, 632.025, 625.798625, 631.953125,
632.155188, 639.860125, 641.00775, 637.589125, 637.24325, 636.01825,
638.107625, 641.851813, 642.569438, 647.145438, 646.820625, 642.397875,
646.79425, 651.9535, 655.33225, 660.207125, 666.556125, 668.464313,
671.055375, 674.749375, 682.068125, 678.8675, 676.1265, 672.481563,
670.703063, 672.062563, 679.603063, 686.332625, 694.170375, 700.403375,
707.086875, 715.055063, 722.802563, 730.246063, 738.165688, 745.08375,
751.146, 748.617188, 746.678813, 746.024563, 742.908813, 745.792625,
750.472063, 754.7605, 756.962, 764.113563, 766.416875, 773.6025,
779.524438, 780.54775, 781.245688, 782.124438, 785.844, 786.494,
793.205188, 798.625563, 793.336563, 788.13475, 783.140375, 779.378563,
777.683813, 780.69875, 787.14325, 779.648813, 772.613438, 769.153125,
769.227688, 765.929875, 763.41925, 758.943313, 758.803438, 755.02475,
758.635313, 755.889125, 749.726063, 741.674063, 735.94675, 729.9425,
729.948, 732.53975, 724.347438, 732.07725, 739.425688, 747.238313,
749.815563, 757.045125, 765.36825, 761.700813, 756.485438, 752.37975,
748.731938, 744.853563, 738.498438, 733.378375, 731.44225, 739.654313,
743.952688, 740.586188, 734.3525, 726.756313, 718.759688, 710.763438,
702.864375, 695.90125, 688.046, 680.367938, 670.204875),
north = c(111.761,
107.051, 105.242, 100.7915, 97.993, 101.1545, 105.5405, 110.7495,
116.397, 120.5795, 119.667, 121.898, 119.9675, 120.177, 115.782,
119.532, 115.8975, 110.499, 117.333, 120.929, 121.4415, 119.9445,
117.39, 110.527, 104.375, 97.7945, 98.896, 97.3895, 94.0785,
94.73, 98.2715, 98.26, 93.482, 90.0065, 83.768, 76.4215, 69.7315,
62.9205, 55.255, 52.0755, 49.4895, 46.7695, 53.727, 59.958, 58.8065,
60.247, 60.1965, 58.7055, 58.782, 61.107, 64.9245, 70.0245, 72.88,
73.6025, 74.997, 75.2305, 77.7445, 78.171, 80.7555, 84.688, 84.1965,
83.564, 83.2995, 82.2695, 80.292, 84.5905, 87.9925, 93.4435,
94.249, 90.7505, 88.96, 92.3015, 96.4505, 103.246, 109.2565,
116.906, 122.0885, 128.354, 136.198, 143.92, 151.001, 158.5445,
163.667, 166.2495, 167.25, 169.1685, 165.457, 173.248, 174.298,
169.808, 168.851, 169.669, 165.8305, 161.0635, 154.883, 161.614,
164.554, 167.945, 176.129, 183.5995, 179.9475, 181.639, 189.343,
190.8625, 188.6405, 192.2815, 196.7395, 201.929, 202.4355, 198.6235,
201.474, 207.0945, 206.12, 207.3115, 205.5975, 210.408, 213.2685,
216.1175, 220.7955, 222.779, 224.8275, 225.788, 228.3695, 228.715,
233.218, 242.7005, 248.682, 249.457, 253.9695, 258.6515, 265.281,
265.8665, 266.8735, 272.053, 279.341, 286.609, 293.7815, 301.3105,
309.0375, 315.479, 323.254, 330.385, 334.1935, 335.4215, 340.129,
345.8245, 353.308, 360.851, 368.452, 375.9125, 383.465, 390.732,
398.7125, 405.3755, 411.065, 415.045, 422.3715, 429.8395, 437.277,
444.9855, 453.0285, 459.9115, 465.8625, 470.151, 474.1165, 477.7425,
480.8905, 485.017, 489.3505, 495.0075, 499.934, 503.1085, 499.2055,
496.7625, 495.3895, 496.865, 500.318, 497.055, 497.705, 494.4745,
496.3685, 491.1985, 491.6935, 489.49, 493.9795, 501.523, 507.9295,
504.9965, 497.152, 503.3365, 502.907, 498.0905, 497.7135, 496.94,
493.1595, 492.287, 489.9165, 492.5375, 493.548, 490.967, 483.0055,
482.201, 490.2815, 491.8335, 484.8545, 482.781, 474.8165, 473.0705,
480.3425, 472.1715, 468.8685, 476.651, 481.32, 474.7985, 475.4645,
479.503, 474.766, 470.4505, 467.3145, 465.864, 463.654, 461.1125,
459.3925, 454.5375, 452.687, 443.8865, 444.3065, 451.4485, 457.7515,
463.096, 470.414, 475.6315, 476.1755, 473.357, 467.075, 459.4815,
459.3755, 464.782, 464.875, 465.694, 462.5555, 464.8, 463.1885,
462.725, 461.2715, 462.521, 465.8655, 464.5585, 460.5085, 456.086,
449.595, 446.604, 444.8285, 437.737, 432.025, 425.496, 422.555,
420.4055, 419.884, 416.7505, 413.4625, 414.2685, 411.1705, 411.021,
403.435, 401.2055, 400.3325, 397.1535, 394.517, 389.9, 382.0995,
375.437, 367.535, 361.227, 355.1195, 349.3295, 342.219, 335.8645,
330.8, 324.334, 317.356, 310.5265, 307.0365, 300.0215, 292.5305,
285.4165, 277.625, 269.8635, 267.0335, 269.963, 269.5365, 273.576,
269.09, 269.743, 269.0615, 267.9115, 267.6705, 267.2245, 261.5795,
255.2725, 247.9015, 240.0395, 232.673, 226.0275, 222.804, 230.096,
237.04, 239.545, 247.0055, 246.4635, 241.737, 235.974, 228.906,
221.1565, 214.5615, 208.1565, 204.8215, 205.266, 199.1955, 193.133,
186.8895, 180.5855, 187.108, 194.386, 199.101, 198.6915, 194.921,
197.2955, 205.304, 200.563, 199.8925, 205.232, 196.9885, 193.748,
188.6695, 181.1865, 182.4225, 182.355, 185.664, 190.334, 185.579,
179.157, 179.706, 177.326, 175.5195, 176.2825, 172.8805, 171.5365,
170.572, 165.084, 159.1015, 152.2355, 145.1195, 138.108, 141.5515,
139.683, 135.9955, 136.1905, 133.044, 125.3875, 126.208, 124.2705,
124.561, 124.8235, 124.122, 121.1025, 122.3335, 121.7645, 111.761
)
), class = "data.frame", row.names = c(NA, -369L)
)
# ---- Variogram ----
geodat_parana = as.geodata(parana, coords.col = 1:2, data.col = 3)
dist_mat = dist(parana[1:2])
min_dist = min(dist_mat)
max_dist = 500
breaks = seq(from = min_dist, to = max_dist, length.out = 50)
emp_variog = variog(geodat_parana, breaks = breaks, option = "bin")
# ---- Matern Variofit ----
# initial values
nugget_0 = 0
sill_0 = 7000
part_sill_0 = sill_0 - nugget_0
range_0 = 400
for(min_kappa in c(0.1, 0.2, 0.3, 0.4, 0.5)){
kappa_list = seq(min_kappa, 50, length.out = 5)
print(paste0("kappa_list = ", kappa_list))
MSE_list = c()
for(kappa_temp in kappa_list){
set.seed(123)
print(paste0("kappa_temp = ", format(kappa_temp, digits = 22, scientific = FALSE)))
kappa_temp <- round(kappa_temp, digits = 1)
variofit_temp = variofit(
emp_variog,
fix.nugget = FALSE, nugget = nugget_0,
fix.kappa = TRUE, kappa = kappa_temp,
ini.cov.pars = c(part_sill_0,range_0),
cov.model = "matern",
messages = F
)
MSE_list = append(variofit_temp$value, MSE_list)
}
print(paste0("MSE for kappa = ", kappa_list[length(kappa_list)],
" is ", MSE_list[length(MSE_list)]))
plot(kappa_list, MSE_list, ylim = c(1.5e+09, 2e+10))
}
Full output :
[1] "kappa_list = 0.1" "kappa_list = 12.575" "kappa_list = 25.05" "kappa_list = 37.525" "kappa_list = 50"
[1] "kappa_temp = 0.1000000000000000055511"
[1] "kappa_temp = 12.57499999999999928946"
[1] "kappa_temp = 25.05000000000000071054"
[1] "kappa_temp = 37.52499999999999857891"
[1] "kappa_temp = 50"
[1] "MSE for kappa = 50 is 21781768624.2714"
[1] "kappa_list = 0.2" "kappa_list = 12.65" "kappa_list = 25.1" "kappa_list = 37.55" "kappa_list = 50"
[1] "kappa_temp = 0.2000000000000000111022"
[1] "kappa_temp = 12.64999999999999857891"
[1] "kappa_temp = 25.09999999999999786837"
[1] "kappa_temp = 37.54999999999999715783"
[1] "kappa_temp = 50"
[1] "MSE for kappa = 50 is 13067343235.6359"
[1] "kappa_list = 0.3" "kappa_list = 12.725" "kappa_list = 25.15" "kappa_list = 37.575" "kappa_list = 50"
[1] "kappa_temp = 0.2999999999999999888978"
[1] "kappa_temp = 12.72500000000000142109"
[1] "kappa_temp = 25.15000000000000213163"
[1] "kappa_temp = 37.57500000000000284217"
[1] "kappa_temp = 50"
[1] "MSE for kappa = 50 is 7316447596.3402"
[1] "kappa_list = 0.4" "kappa_list = 12.8" "kappa_list = 25.2" "kappa_list = 37.6" "kappa_list = 50"
[1] "kappa_temp = 0.4000000000000000222045"
[1] "kappa_temp = 12.80000000000000071054"
[1] "kappa_temp = 25.19999999999999928946"
[1] "kappa_temp = 37.60000000000000142109"
[1] "kappa_temp = 50"
[1] "MSE for kappa = 50 is 4028582940.74542"
[1] "kappa_list = 0.5" "kappa_list = 12.875" "kappa_list = 25.25" "kappa_list = 37.625" "kappa_list = 50"
[1] "kappa_temp = 0.5"
[1] "kappa_temp = 12.875"
[1] "kappa_temp = 25.25"
[1] "kappa_temp = 37.625"
[1] "kappa_temp = 50"
[1] "MSE for kappa = 50 is 2728655320.43029"
The problem is with your reporting. When you print MSE for kappa = 50 is 21781768624.271 you're actually printing the MSE value for some other kappa value, because you're putting new values at the start of the list, not the end of the list, in this expression:
MSE_list = append(variofit_temp$value, MSE_list)
To fix things, just use
MSE_list = c(MSE_list, variofit_temp$value)
(or use append() in place of c(), they do the same thing here.)