A general program in the book :
André Dauphiné, Geographical models with Mathematica
country = "France";
ny = ToExpression@
DialogInput[
DynamicModule[{name = ""},
Column[{"Combien de villes ?", InputField[Dynamic[name], String],
ChoiceButtons[{DialogReturn[name], DialogReturn[]}]}]]];
coord = Take[
Table[Reverse[CityData[c, "Coordinates"]], {c,
CityData[{All, country}]}], ny];
coordxy =
GeoGridPosition[
GeoPosition[#, "WGS84"], {"UTMZone31",
"CentralScaleFactor" -> 0.9996,
"GridOrigin" -> {500000, 0}}][[1]] & /@ coord;
pop = Take[
Table[QuantityMagnitude[CityData[c, "Population"]], {c,
CityData[{All, country}]}], ny];
nom = Take[Table[CityData[c, "Name"], {c, CityData[{All, country}]}],
ny];
n1 = Length[coord];
xyz = Partition[Flatten[Riffle[coordxy, pop]], 3];
minx = Min[xyz[[All, 1]]];
miny = Min[xyz[[All, 2]]];
maxx = Max[xyz[[All, 1]]];
maxy = Max[xyz[[All, 2]]];
Print["Paramètres du modèle d'ordre 3"]
ln = GeneralizedLinearModelFit[
xyz, {x, y, x*y, x^2, y^2, x*y^2, y*x^2, x^3, y^3}, {x, y}];
Normal[ln]
test = ln[{"AIC", "BIC"}];
Grid[{{test_AIC, test[[1]]}, {test_BIC, test[[2]]}}, Frame -> All]
fit = ln["BestFit"];
resid = ln[{"FitResiduals", "StandardizedPearsonResiduals",
"CookDistances"}];
xyzr = Partition[Flatten[Riffle[nom, Transpose[resid]]], 4];
nn = {Noms_ville, Résidus_bruts, Résidus_standardisés,
Distance_Cook};
xyzr = Insert[xyzr, nn, 1];
Print[]
Print["Surface de tendance et valeurs des données"]
surf = Show[ContourPlot[fit, {x, minx, maxx}, {y, miny, maxy}],
Frame -> False, ClippingStyle -> Automatic,
ColorFunction -> ColorData[{"GrayLevel", "Reverse"}]]
GeoGraphics[{GeoStyling[{"GeoImage",
ImageCrop[Rasterize[surf, RasterSize -> 400]]}], EdgeForm[Thin],
Polygon[Entity["Country", country]]},
GeoRange -> {{40.5, 52.5}, {-5, 9.5}}]
Print[]
Grid[xyzr, Alignment -> ".", Frame -> All,
Background -> {{Yellow, LightGray, LightBlue, LightRed}, None}]