Nick, thanks for sharing this with us. The good news of course is that
both gstat and geoR give the same results for isotropic kriging. The trouble is that software developers like Paulo and me are mainly bothered with developing our own thing, and not too much with compatibility to the others, although last March we agreed to make an effort for exactly this reason. A similar thing that came up a few years ago on ai-geostats was the specification of the parameters for the exponential and Gaussian variogram models -- see http://ai-geostats.org/Geostats_Faq/FAQ_software_conventions.htm and we need such a page for modelling (and fitting!) anisotropy parameters. I did a little experiment with two generated anisotropic random fields, x and y, R file and png images attached. The geoR output does estimate the ratio fairly well: a value of 0.3 for gstat corresponds to 3.333.. in geoR; however I do not understand the psiA value; geoR docs say: Anisotropy angle: defined here as the azimuth angle of the direction with greater spatial continuity, i.e. the angle between the y-axis and the direction with the maximum range. another difference is that gstat uses degrees (0-180), whereas I think geoR uses radial degrees (?), (0-pi). Still, in that case geoR should give values close to 0.11 (20/180) and 0.611 (110/180). Paulo? Also, the sigmasq values worry me a bit: > var(x[,3]) [1] 0.8750844 > var(y[,3]) [1] 0.9214435 generating variogram(x): 1sph(10, 20, 0.3); > likfit(as.geodata(x),ini=c(1,10),fix.nugget=T,cov.model="spherical",fix.psiA=F, psiA=0.111, fix.psiR=F, psiR=3) .... likfit: estimated model parameters: beta sigmasq phi psiA psiR "-0.0713" " 2.5342" " 9.5813" " 0.3248" " 3.7673" variogram(x): 1sph(10, 110, 0.3); > likfit(as.geodata(y),ini=c(1,10),fix.nugget=T,cov.model="spherical",fix.psiA=F, psiA=0.05*c(1:20)*pi, fix.psiR=F, psiR=3) likfit: searching for best initial value ... selected values: sigmasq phi tausq kappa lambda psiR psiA initial.value " 1.0" "10.0" " 0.0" " 0.5" " 1.0" " 3.0" " 2.0" status "est" "est" "fix" "fix" "fix" "est" "est" likelihood value: -539.58582544843 .... beta sigmasq phi psiA psiR " 0.1186" " 2.3785" "10.0816" " 1.9426" " 3.1296" likfit: maximised log-likelihood = -337.4 I hope we can resolve this soon. As a short-term answer: you can obtain block mean values by point kriging on a fine grid, and then averaging point kriged values within a block. It does not give you block kriging standard errors, though. -- Edzer Nick Hamm wrote: Hi I have a set of two dimensional data( measured on a grid). Initial analysis strongly suggests that there is directional anisotropy, so I've been trying to deal with this.One way I've found is to fit the variogram model by maximum likelihood (or REML) using geoR. geoR seems to account for anisotropy by deforming the locations, according to the anisotropy parameters (direction of maximum range, and the anisotropy ratio). Cross--validation suggests that this performs well, relative to fitting various models to isotropic sample variogram. The thing is, I want to use gstat for kriging. The reason for this is that I want to do block kriging. However, gstat seems to account for anisotropy by adjusting the variogram range according to the anisotropy parameters. So it seems that geoR deforms the coordinates, whereas gstat "deforms" the variogram.... The upshot of this is that, when I use gstat to do punctual kriging using the parameters estimated in geoR, I get different results to what I get when I use geoR for kriging. A possible alternative, is to use geoR to predict the deformed locations, and then use gstat to do the kriging WITHOUT specifying the anisotropy parameter. If I do this then gstat and geoR give the same results for kriging. This suggests that there is no difference in the algorithms used for kriging between the two packages. This approach seems to be OK for doing punctual kriging... However, I want to do block kriging (square blocks). The approach I've adopted for punctual kriging seems to fall over here. This is because I would specify the block size (in gstat) in the deformed coordinate space, whereas I actually want the blocks to be specified in the original geographic coordinate space. Am I making sense? What a mess!!! Can any one suggest an alternative, or is this all a bit silly? cheers Nick |
<<inline: x.png>>
<<inline: y.png>>
"x" <- structure(list(x = c(-9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5), y = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5), attr = c(-0.976396, 0.699642, 1.02706, 1.8274, -0.0842941, 0.730785, 1.17241, 0.17138, -0.499801, 0.0290917, 0.971627, 0.159505, 0.0978222, -0.0628971, 1.10981, 1.60993, 0.157566, -1.31656, -0.814423, -1.08538, -0.664614, 1.3288, 1.7715, 0.648596, 0.512488, 0.125008, 1.15626, -0.293581, -1.40186, 0.345669, 1.36306, 0.828358, -0.342381, 0.00905169, 1.24502, 0.576963, -0.0811922, -1.44835, -0.319535, -0.341521, 0.515442, 1.15643, 1.48897, 0.418224, 0.129235, 0.0275096, -0.505079, -1.1466, -0.922303, 0.501884, 1.3434, 0.690577, -0.0654742, 0.240356, 0.833001, 0.465874, -1.25548, -0.841794, 0.228543, 0.924724, 0.461841, 2.01067, 0.684979, -0.194078, -0.311548, 0.73425, -0.682227, -1.32422, -1.31838, 0.762056, 1.92055, 0.664635, 0.169215, 1.13904, -0.0118303, -0.42341, -2.09248, -0.443041, -0.521656, 1.20924, 1.60895, 0.541853, -1.32658, -0.34973, 0.320017, 0.0198023, -1.85797, -0.837399, -0.861881, 0.752807, 1.19442, 0.55206, -0.337885, -0.268511, -0.0354799, -0.577635, -1.1632, -0.0418927, -0.190198, 1.38387, 0.864981, -0.230297, -0.908201, -0.317284, 0.540593, 0.588873, -2.15811, -1.29086, 0.564227, 0.45145, 0.00254214, 0.223103, -0.72224, -0.691557, 0.293142, -0.229129, -0.0911508, -0.446708, 0.0914658, 1.91125, -0.0655018, -0.409656, -0.785249, 0.0342057, 1.0731, 0.226422, -1.59694, -1.20873, 0.506341, 0.979477, 0.164993, -0.676357, -0.608324, -0.989886, -0.74416, 0.844582, -0.828441, 0.627724, 1.75291, 1.97964, -0.147929, -0.929661, -1.73665, -0.349646, 1.48683, -0.289314, -2.0301, -0.418626, 0.404249, 1.43853, 0.0772116, -1.2291, -2.07862, -1.02916, -0.3847, 0.233113, -0.484702, 1.29267, 2.53204, 1.27968, -0.598517, -0.369162, -0.592886, 0.937435, 1.3464, -0.48181, -0.404962, 0.930906, 1.7464, 1.18944, -0.177596, -1.65704, -1.86563, -0.59391, 0.396574, 0.285294, 0.150976, 1.0885, 3.13262, 1.42033, -1.07894, -0.422835, 0.340313, 1.04254, 1.02172, -0.930621, -0.130838, 0.891712, 1.83336, 0.903686, -0.121988, -1.51042, -0.680568, 0.262939, 0.345589, -0.0345113, 1.29514, 1.91254, 1.89058, 1.40523, -0.352911, -0.471958, 1.17492, 1.06129, 0.0261414, -0.00657254, 0.270249, 0.831697, 1.82552, 0.873856, -0.684648, -1.86114, -0.119203, 0.311865, 0.167739, -0.0225983, 1.37158, 1.93348, 0.282233, 1.81257, -0.0534109, -0.178551, 0.559363, 0.706015, -0.252069, -0.175302, -0.720996, -0.166321, 0.743316, 0.812269, -1.17577, -0.932122, -0.566561, 0.718219, 0.188124, 0.735218, 1.6418, 0.178471, 0.407502, 1.36858, -0.000146809, 0.291759, -0.468605, 0.0206043, -0.944232, -0.61307, 0.442481, -0.968002, 1.81035, 0.181125, -1.23142, -0.763612, 0.0737858, 1.12902, -0.254041, 0.927865, 0.886834, -0.571706, -0.13527, 1.27622, -0.907889, -0.657191, -0.758956, 0.561848, -0.681267, -0.222691, 0.463942, 0.224039, 1.25724, 0.25348, -1.35065, -0.709735, 1.73603, 0.84702, -0.958807, -0.203773, -0.749048, -0.215812, 0.948336, 0.8805, -0.953084, -1.49239, -0.0665188, -0.533943, -0.0828005, -0.870736, 0.65511, -0.137514, 0.552696, 0.0139825, -1.84643, -0.360674, 2.06584, 0.0382007, -0.198137, -0.187399, -1.5025, -0.560836, 0.790284, 0.477532, -0.86596, -0.68644, -0.690822, 0.627588, -0.635185, 0.685255, 0.0380479, 0.681636, 0.503216, -0.0772367, -0.858596, 0.890215, 0.857633, -0.705095, -0.759286, -1.27039, -0.652695, -0.123117, 1.34745, -0.0751269, -0.418802, -0.0673295, -0.265125, 1.11624, -1.15422, 0.19013, -0.219525, 0.53884, -0.0382443, -0.948078, -0.477328, -0.634133, -0.369468, -1.50007, -1.08964, -1.31225, 0.270723, 0.442437, 1.9753, 0.434132, -1.16969, -0.198799, 0.606902, 0.711261, -0.719003, 0.473304, -0.649792, 0.472475, -0.314922, -1.41962, -0.428929, -0.653978, -0.777673, -1.13312, -1.23852, -0.990636, 0.0964782, 0.26112, 1.37736, 0.977292, 0.0433025, 0.878182, 1.06444, 0.89754, -0.130647, -0.538042, 0.0716454, 0.191813, -0.884949, -2.00202, -0.0712418, -1.35596, -0.373052, -1.25609, -1.68115, -1.31557, -0.605708, 0.638957, 0.841501, 0.701605, -0.322548, -0.374385, 0.692801, 0.963631, -0.0116963, -0.122099, 0.125823, -1.04663, -0.769083, -1.01459, -0.524525, -1.90085, -0.898233, -1.45772, -1.68847, -1.82645, -0.429825, 0.0917008, 0.948507, 0.305031)), .Names = c("x", "y", "attr"), class = "data.frame", row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "120", "121", "122", "123", "124", "125", "126", "127", "128", "129", "130", "131", "132", "133", "134", "135", "136", "137", "138", "139", "140", "141", "142", "143", "144", "145", "146", "147", "148", "149", "150", "151", "152", "153", "154", "155", "156", "157", "158", "159", "160", "161", "162", "163", "164", "165", "166", "167", "168", "169", "170", "171", "172", "173", "174", "175", "176", "177", "178", "179", "180", "181", "182", "183", "184", "185", "186", "187", "188", "189", "190", "191", "192", "193", "194", "195", "196", "197", "198", "199", "200", "201", "202", "203", "204", "205", "206", "207", "208", "209", "210", "211", "212", "213", "214", "215", "216", "217", "218", "219", "220", "221", "222", "223", "224", "225", "226", "227", "228", "229", "230", "231", "232", "233", "234", "235", "236", "237", "238", "239", "240", "241", "242", "243", "244", "245", "246", "247", "248", "249", "250", "251", "252", "253", "254", "255", "256", "257", "258", "259", "260", "261", "262", "263", "264", "265", "266", "267", "268", "269", "270", "271", "272", "273", "274", "275", "276", "277", "278", "279", "280", "281", "282", "283", "284", "285", "286", "287", "288", "289", "290", "291", "292", "293", "294", "295", "296", "297", "298", "299", "300", "301", "302", "303", "304", "305", "306", "307", "308", "309", "310", "311", "312", "313", "314", "315", "316", "317", "318", "319", "320", "321", "322", "323", "324", "325", "326", "327", "328", "329", "330", "331", "332", "333", "334", "335", "336", "337", "338", "339", "340", "341", "342", "343", "344", "345", "346", "347", "348", "349", "350", "351", "352", "353", "354", "355", "356", "357", "358", "359", "360", "361", "362", "363", "364", "365", "366", "367", "368", "369", "370", "371", "372", "373", "374", "375", "376", "377", "378", "379", "380", "381", "382", "383", "384", "385", "386", "387", "388", "389", "390", "391", "392", "393", "394", "395", "396", "397", "398", "399", "400")) "y" <- structure(list(V1 = c(-9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5), V2 = c(9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 9.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 8.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 7.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 6.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 4.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 3.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 2.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 1.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -0.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -1.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -2.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -3.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -4.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -5.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -6.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -7.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -8.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5, -9.5), V3 = c(-0.358227, -0.883175, 0.428496, -0.213165, -0.121203, 0.504585, 1.29542, 0.644304, 0.350139, 0.114496, 0.174196, 0.328178, -0.0868328, 0.226939, 0.346281, 0.401713, 1.30887, 1.76147, 0.853352, 0.121347, -0.337769, -0.953408, -1.69177, -1.00548, -0.837798, -0.422826, -0.0688711, 0.323329, 0.995217, 0.638197, 0.762255, 0.870313, 0.711054, 1.04473, 1.34499, 1.17326, 0.218925, 0.936686, 0.861577, 0.383784, 0.938214, 1.08439, -0.0492151, -0.251227, -1.49285, -1.03058, -2.03424, -1.24326, -0.634366, -0.193798, 0.390262, 1.66108, 1.92824, 1.5671, 0.191612, 0.338176, 1.50281, 1.0535, 0.551203, 0.0413718, 0.377746, 0.723457, 0.142913, -0.89396, -0.782161, -0.981936, -0.973828, -1.30591, -1.73934, -1.84269, -0.0832526, 1.38839, 1.51343, 2.0416, 2.24402, 1.90511, 1.97194, 0.324144, -0.23269, -0.692648, 1.05821, 1.03786, 0.784226, -0.339956, -0.676398, -1.01696, -0.95513, -0.487205, -1.22426, -1.56292, -1.04767, -0.265129, 0.585775, 1.53192, 1.06778, 2.4384, 1.06568, 1.29809, 0.924595, -0.606835, 1.10346, 1.56817, 1.02077, 1.24847, 0.248581, 0.224182, -0.717934, -0.188762, -0.3741, -1.05474, -0.597772, -0.564637, -0.921164, -0.271499, -0.906596, 0.492014, 0.630646, 0.854272, 1.96541, 0.999921, 0.496202, 0.351973, 1.60986, 2.33434, 1.51661, 0.167276, -0.146389, -0.225696, -0.892759, -0.0899573, -1.39838, -0.981012, -0.58096, -0.694041, -1.8745, -1.89785, -1.75994, -1.4497, 0.201812, 0.619787, -0.485159, 0.273018, 1.14061, 1.97519, 2.11462, 0.839004, 0.717381, 0.824517, 0.259545, 0.817707, -1.27184, -0.457396, 0.338396, -0.228584, 0.194588, -0.077346, -0.926435, -1.50319, -1.55432, -1.17051, -1.09185, -0.917399, -0.23146, -0.375039, 0.155894, 0.220507, 0.231177, 0.97191, -0.0836007, -1.33339, -0.622234, -0.201292, -0.854076, -0.367267, 0.384303, 1.13566, 1.03534, 0.516219, -0.48998, -1.18473, 0.196995, 0.262778, -0.0359566, 0.47255, -0.0905125, 0.357794, -0.285613, -0.00798318, -0.11457, -0.761042, -0.439379, -0.835278, -0.639411, -0.785074, -0.00523192, 0.1666, 1.05992, 1.06818, 0.833633, 0.385857, 0.921412, 0.276106, 0.463932, 0.358063, -0.0233177, -0.238753, -0.42071, -0.251407, -0.282355, -0.0486227, -0.130755, -0.595462, 0.490839, 0.339338, -0.0560451, 0.489945, 0.517167, 0.40683, 1.00145, 0.589293, 1.53035, 0.337886, 0.536092, 0.748265, 0.985178, 0.618256, -0.744364, -1.40118, -0.225276, -0.133966, -0.741826, -0.662154, -1.13363, -0.199749, 0.218788, -0.381617, 0.579232, 1.01703, 0.765987, 0.745737, 1.26144, 0.722743, 0.232609, 0.296584, 0.867954, 2.10431, 1.52535, 1.06925, 0.314842, 0.356067, 0.0942154, -0.108979, -0.988771, -1.90368, -1.91997, -1.34157, -0.217298, 1.22345, 1.15753, 0.441481, 1.55639, 0.882917, -0.792265, -0.0559937, 0.201641, 1.13765, 1.73539, 1.28896, 0.854913, 0.901594, -0.0555984, -0.7484, -0.195218, -0.602409, -1.71683, -2.07703, -1.95228, -0.548744, 0.0371549, 0.595531, 1.21788, 1.17501, -0.147261, -0.510002, -1.14883, -0.700317, -0.134978, 0.268485, 0.327767, 0.1389, 0.243102, -0.289124, 0.309648, -0.290261, -0.564847, -0.641719, -1.75654, -2.05613, -1.75529, -0.0565981, 1.88894, 1.48677, 0.646488, 0.0356698, -0.113009, -0.45765, -0.720445, -0.786172, -0.724996, 0.0962278, 0.0526066, -2.20774, -1.21379, -0.291403, -0.665068, -0.912585, -0.64879, -1.8832, -1.095, -0.90682, 0.646028, 2.11661, 2.45733, 1.26387, 1.07263, -0.416093, -0.0959732, -0.327427, -0.167697, 0.422089, 0.144705, -0.401572, -0.348742, -0.650393, -0.136084, -1.01573, -0.546468, -1.49028, -0.893513, -0.51415, 0.0848532, 1.01496, 2.22971, 2.30453, 2.29134, 0.976251, 0.176245, -0.561158, -0.930977, -0.628719, -0.172885, 0.518523, -0.825708, -1.04278, -0.0333577, -1.53102, -0.990126, -0.942557, -1.55075, -0.831715, 0.198945, 0.367643, 1.02887, 1.86658, 0.936428, 1.05558, 0.426798, 0.176225, 0.113876, -0.509315, -0.00407159, -0.0213874, 0.444478, 0.545452, -0.359192, -0.562048, -0.0992335, -0.869162, -1.66234, -1.9709, -0.449777, -0.263866, -0.270827, 0.194169, 1.67132, 1.44206, 1.94036, 0.553896, -0.106475, 0.304114, 0.0289134, 0.0534254, -0.0438635, -0.445235, -0.455341, 0.202635, -0.232922, -0.798149, 0.756377, -0.150466)), .Names = c("V1", "V2", "V3"), class = "data.frame", row.names = c("1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", "21", "22", "23", "24", "25", "26", "27", "28", "29", "30", "31", "32", "33", "34", "35", "36", "37", "38", "39", "40", "41", "42", "43", "44", "45", "46", "47", "48", "49", "50", "51", "52", "53", "54", "55", "56", "57", "58", "59", "60", "61", "62", "63", "64", "65", "66", "67", "68", "69", "70", "71", "72", "73", "74", "75", "76", "77", "78", "79", "80", "81", "82", "83", "84", "85", "86", "87", "88", "89", "90", "91", "92", "93", "94", "95", "96", "97", "98", "99", "100", "101", "102", "103", "104", "105", "106", "107", "108", "109", "110", "111", "112", "113", "114", "115", "116", "117", "118", "119", "120", "121", "122", "123", "124", "125", "126", "127", "128", "129", "130", "131", "132", "133", "134", "135", "136", "137", "138", "139", "140", "141", "142", "143", "144", "145", "146", "147", "148", "149", "150", "151", "152", "153", "154", "155", "156", "157", "158", "159", "160", "161", "162", "163", "164", "165", "166", "167", "168", "169", "170", "171", "172", "173", "174", "175", "176", "177", "178", "179", "180", "181", "182", "183", "184", "185", "186", "187", "188", "189", "190", "191", "192", "193", "194", "195", "196", "197", "198", "199", "200", "201", "202", "203", "204", "205", "206", "207", "208", "209", "210", "211", "212", "213", "214", "215", "216", "217", "218", "219", "220", "221", "222", "223", "224", "225", "226", "227", "228", "229", "230", "231", "232", "233", "234", "235", "236", "237", "238", "239", "240", "241", "242", "243", "244", "245", "246", "247", "248", "249", "250", "251", "252", "253", "254", "255", "256", "257", "258", "259", "260", "261", "262", "263", "264", "265", "266", "267", "268", "269", "270", "271", "272", "273", "274", "275", "276", "277", "278", "279", "280", "281", "282", "283", "284", "285", "286", "287", "288", "289", "290", "291", "292", "293", "294", "295", "296", "297", "298", "299", "300", "301", "302", "303", "304", "305", "306", "307", "308", "309", "310", "311", "312", "313", "314", "315", "316", "317", "318", "319", "320", "321", "322", "323", "324", "325", "326", "327", "328", "329", "330", "331", "332", "333", "334", "335", "336", "337", "338", "339", "340", "341", "342", "343", "344", "345", "346", "347", "348", "349", "350", "351", "352", "353", "354", "355", "356", "357", "358", "359", "360", "361", "362", "363", "364", "365", "366", "367", "368", "369", "370", "371", "372", "373", "374", "375", "376", "377", "378", "379", "380", "381", "382", "383", "384", "385", "386", "387", "388", "389", "390", "391", "392", "393", "394", "395", "396", "397", "398", "399", "400"))