require(geoR) data(parana) plot(parana) ml <- likfit(parana, ini=c(100, 40), trend="1st") gr <- pred_grid(parana$borders, by=10) points(parana) points(gr, pch=19, cex=0.3, col=2) args(krige.conv) args(output.control) kc <- krige.conv(parana, loc=gr, krige=krige.control(obj.model=ml, trend.l="1st"), output=output.control(simulations.p=1000)) names(kc) image(kc) image(kc, value=kc$krige.var) dim(kc$simulations) COL <- gray(seq(1, 0.2, length=21)) par(mfrow=c(2,3), mar=c(3,3,1,1)) image(kc, value=kc$simulations[,1], col=COL) image(kc, value=kc$simulations[,2], col=COL) image(kc, value=kc$simulations[,3], col=COL) #image(kc, value=kc$simulations[,4], col=COL) image(kc, col=COL) image(kc, value=apply(kc$simulations, 1, mean), col=COL) plot(kc$pred, apply(kc$simulations, 1, mean), asp=1); abline(0,1) summary(kc$pred - apply(kc$simulations, 1, mean)) p400 <- apply(kc$simulations, 1, function(x) mean(x > 400)) p300 <- apply(kc$simulations, 1, function(x) mean(x > 300)) p200 <- apply(kc$simulations, 1, function(x) mean(x > 200)) length(p300) par(mfrow=c(2,3), mar=c(3,3,1,1)) image(kc, value=p300, x.leg=c(200, 700), y.leg=c(-50, 0)) image(kc, value=p200, x.leg=c(200, 700), y.leg=c(-50, 0)) image(kc, value=p400, x.leg=c(200, 700), y.leg=c(-50, 0)) image(kc, value=p300, x.leg=c(200, 700), y.leg=c(-50, 0), breaks=c(0, 0.2, 0.4, 0.8, 1), col=heat.colors(4)) image(kc, value=p300, x.leg=c(200, 700), y.leg=c(-50, 0), breaks=c(0, 0.2, 0.4, 0.8, 1), col=gray(c(1, 0.3, 0.6, 0))) prop300 <- apply(kc$simulations, 2, function(x) mean(x > 300)) length(prop300) hist(prop300, prob=T);lines(density(prop300));rug(prop300) ## v <- variog(parana, trend="1st", max.dist=450 v10 <- variog(parana, trend="1st", uvec=seq(0, 450, length=10),max.dist=450) plot(v10) ## ef <- eyefit(v) vf <- variofit(v) lines(vf, col=2) vf1 <- variofit(v, cov.model="circ") lines(vf1, col=2)