###################################### ### Cviceni 19. 3. 2014 - nahodna pole ###################################### ############### # Isinguv model ############### # java applet s MCMC simulaci: # http://physics.ucsc.edu/~peter/ising/ising.html # # Init cold, Init hot, Init warm - pocatecni konfigurace simulace # Magnetization - prumer orientaci vsech spinu ################ # knihovna spdep ################ library(spdep) data(package="spdep") # ruzne priklady regionalnich dat library(maptools) example(columbus) # data z 49 oblasti v Columbusu (Ohio) plot(columbus) coords <- coordinates(columbus) points(coords,pch=19) ### sousede col.gal.nb # relace sousedstvi plot(col.gal.nb,coords,add=TRUE) summary(col.gal.nb) is.symmetric.nb(col.gal.nb) print(cards <- card(col.gal.nb)) # pocty sousedu ### vahy col.w <- nb2listw(col.gal.nb,style="W") # vahy normalizovane na soucet 1 col.w$weights col.b <- nb2listw(col.gal.nb,style="B") # binarni vahy col.b$weights ### druhy stupen sousedstvi col.lags <- nblag(col.gal.nb, 2) # sousede sousedu plot(col.lags[[2]], coords, add=TRUE, col="red", lty=2) summary(col.lags[[2]], coords) ### k nejblizsich sousedu library(RANN) col.k4.nb <- knn2nb(knearneigh(coords, 4)) # 4 nejblizsi is.symmetric.nb(col.k4.nb) plot(columbus,border="grey") plot(col.k4.nb,coords,arrows=TRUE,add=TRUE) ### rozdil mezi dvema relacemi sousedstvi col.k1.nb <- knn2nb(knearneigh(coords, 1)) col.k2.nb <- knn2nb(knearneigh(coords, 2)) diffs <- diffnb(col.k2.nb,col.k1.nb) plot(columbus, border="grey") plot(col.k1.nb,coords,add=TRUE) plot(diffs,coords,add=TRUE,col="red",lty=2) ### sousedstvi definovane na zaklade euklidovske vzdalenosti col.d.nb <- dnearneigh(coords,0,0.5) # sousede, pokud jsou ve vzdalenosti mezi 0 a 0.5 is.symmetric.nb(col.d.nb) plot(columbus,border="grey") plot(col.d.nb,coords,add=TRUE) ### Delaunayova triangulace col.tri.nb <- tri2nb(coords) plot(columbus,border="grey") plot(col.tri.nb,coords,add=TRUE) # zadny vrchol nelezi uvnitr nejake kruznice opsane trojuhelniku ### spojita data columbus$CRIME # vloupani a kradeze aut na 1000 domacnosti moran.test(columbus$CRIME,listw=col.w) # Moranuv index pro normalizovane vahy, test za predp. znahodneni moran.test(columbus$CRIME,listw=col.w,randomisation=FALSE) # Moranuv index pro normalizovane vahy, test za predp. normality moran.test(columbus$CRIME,listw=col.b) # Moranuv index pro binarni vahy moran.test(columbus$CRIME,listw=nb2listw(col.k4.nb)) # Moranuv index pro normalizovane vahy dane 4 nejblizsimi sousedy geary.test(columbus$CRIME,listw=col.w) # Gearyho index pro normalizovane vahy, test za predp. znahodneni geary.test(columbus$CRIME,listw=col.b,randomisation=FALSE) # Gearyho index pro normalizovane vahy, test za pred. normality geary.test(columbus$CRIME,listw=col.b) # Gearyho index pro binarni vahy geary.test(columbus$CRIME,listw=nb2listw(col.k4.nb)) # Gearyho index pro normalizovane vahy dane 4 nejblizsimi sousedy ### binarni data HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high")) plot(columbus) points(coords[HICRIME=="high",],col="red",pch=19) points(coords[HICRIME=="low",],col="blue",pch=19) joincount.test(HICRIME, listw=col.w) # BB statistika joincount.test(HICRIME, listw=col.b) joincount.test(HICRIME, listw=nb2listw(col.k4.nb)) ### sousede pro pravidelnou mrizku nb7rt <- cell2nb(7, 7) nb7rt summary(nb7rt) xyc <- attr(nb7rt, "region.id") xy <- matrix(as.integer(unlist(strsplit(xyc, ":"))), ncol=2, byrow=TRUE) plot(nb7rt, xy) print(card(nb7rt)) nb7rtt <- cell2nb(7, 7, torus=TRUE) summary(nb7rtt) print(card(nb7rtt)) ### nekorelovana data data7rt <- runif(49,0,1) moran.test(data7rt,nb2listw(nb7rt)) geary.test(data7rt,nb2listw(nb7rt)) ####################### # knihovna RandomFields ####################### library(RandomFields) ### R >= 3.0.2 ### simulace gaussovskych nahodnych poli model1 <- RMspheric(var=2,scale=0.2) # sfericka autokovariancni funkce plot(model1,xlim=c(-0.5,0.5),ylim=c(0,2)) x <- seq(0,1,by=0.02); y <- seq(0,1,by=0.02) # dvourozmerna pravidelna mriz bodu v [0,1]^2 Z1 <- RFsimulate(model1,x,y,grid=TRUE) # realizace stacionarniho gaussovskeho nahodneho pole plot(Z1) z <- seq(0,1,by=0.02) model2 <- RMexp(var=1,scale=0.1)+RMnugget(var=0.2)+RMtrend(mean=0.5) # do modelu lze pripad zbytkovy rozptyl a stredni hodnotu # nyni pouzivame exponencialni autokovariancni funkci plot(model2,ylim=c(0,1.2)) plot(model2,ylim=c(0,1.2), fct.type="Variogram") # semivariogram plot(model2,xlim=c(-0.3,0.3), fct.type="Variogram", dim=2) # semivariogram 2D Z2 <- RFsimulate(model2,x,y,z,grid=TRUE) # simulace 3D gaussovskeho nahodneho pole plot(Z2,MARGIN=c(2,3),MARGIN.slices=1,n.slices=4,col=terrain.colors(64)) # Ukol: vyzkousejte ruznou volbu parametru a pozorujte, jak ovlivnuje vzhled vyslednych realizaci model3 <- RMfbm(alpha=1) # frakcionalni Brownuv pohyb, Hurstuv parametr H=alpha/2 plot(model3,ylim=c(0,1)) # semivariogram - neni omezena funkce Z3 <- RFsimulate(model3,x=seq(0,10,by=0.02),grid=TRUE) plot(Z3) model4 <- RMmatern(nu=2) # cim vetsi nu, tim hladsi pole (nu=0.5 odpovida exponencialni autokovariancni funkci) Z4 <- RFsimulate(model4,x=cbind(x,y),grid=TRUE) plot(Z4) ### empiricky semivariogram ev4 <- RFempiricalvariogram(data=Z4) plot(ev4) plot(ev4, model=model4) ############### # knihovna geoR ############### library(geoR) # knihovna pro geostatistickou analyzu ### Kattegat data(kattegat) # mereni slanosti v uzine Kattegat plot(c(550,770),c(6150,6420),type="n",xlab="X UTM",ylab="Y UTM") points(kattegat, pch=19, add=TRUE) lapply(kattegat$dk, lines, lwd=2) ### Maza data(meuse.riv) # data z knihovny sp plot(meuse.riv,type="l",asp=1) # reka Maza data(meuse) # koncentrace tezkych kovu v pude points(meuse) summary(meuse) plot(y~x,meuse) plot(meuse[,1:4]) # souradnice, kadmium a med meuse_zn <- as.geodata(meuse,data.col=6) # zinek plot(c(178000,182000),c(33000,33400),type="n") points(meuse_zn, pch=19) ### Odhad semivariogramu odhad_gamma <- variog(meuse_zn,max.dist=1500) # neparametricky odhad semivariogram pro zinek plot(odhad_gamma,type="b") smooth <- variog(meuse_zn,option="smooth",max.dist=1500,n.points=100,kernel="normal",band=250) lines(smooth,type="l",lty=2) # jadrovy odhad eyefit(odhad_gamma) # prolozeni parametrickeho modelu ### Neizotropie odhad_gamma_45 <- variog(meuse_zn,max.dist=1500,direction=pi/4) plot(odhad_gamma_45,type="b")