R : Copyright 2005, The R Foundation for Statistical Computing Version 2.1.1 (2005-06-20), ISBN 3-900051-07-0 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for a HTML browser interface to help. Type 'q()' to quit R. > ### *
> ### > attach(NULL, name = "CheckExEnv") > assign(".CheckExEnv", as.environment(2), pos = length(search())) # base > ## add some hooks to label plot pages for base and grid graphics > setHook("plot.new", ".newplot.hook") > setHook("persp", ".newplot.hook") > setHook("grid.newpage", ".gridplot.hook") > > assign("cleanEx", + function(env = .GlobalEnv) { + rm(list = ls(envir = env, all.names = TRUE), envir = env) + RNGkind("default", "default") + set.seed(1) + options(warn = 1) + delayedAssign("T", stop("T used instead of TRUE"), + assign.env = .CheckExEnv) + delayedAssign("F", stop("F used instead of FALSE"), + assign.env = .CheckExEnv) + sch <- search() + newitems <- sch[! sch %in% .oldSearch] + for(item in rev(newitems)) + eval(substitute(detach(item), list(item=item))) + missitems <- .oldSearch[! .oldSearch %in% sch] + if(length(missitems)) + warning("items ", paste(missitems, collapse=", "), + " have been removed from the search path") + }, + env = .CheckExEnv) > assign("..nameEx", "__{must remake R-ex/*.R}__", env = .CheckExEnv) # for now > assign("ptime", proc.time(), env = .CheckExEnv) > grDevices::postscript("RandomFields-Examples.ps") > assign("par.postscript", graphics::par(no.readonly = TRUE), env = .CheckExEnv) > options(contrasts = c(unordered = "contr.treatment", ordered = "contr.poly")) > options(warn = 1) > library('RandomFields') > > assign(".oldSearch", search(), env = .CheckExEnv) > assign(".oldNS", loadedNamespaces(), env = .CheckExEnv) > cleanEx(); ..nameEx <- "CondSimu" > > ### * CondSimu > > flush(stderr()); flush(stdout()) > > ### Name: CondSimu > ### Title: Conditional Simulation > ### Aliases: CondSimu > ### Keywords: spatial > > ### ** Examples > > ## creating random variables first > ## here, a grid is chosen, but any arbitrary points for which > ## data are given are fine. Indeed if the data are given on a > ## grid, the grid has to be expanded before calling `CondSimu', > ## see below. > ## However, locations where values are to be simulated, > ## should be given in form of a grid definition whenever > ## possible > param <- c(0, 1, 0, 1) > model <- "exponential" > > RFparameters(PracticalRange=FALSE) > p <- 1:7 > data <- GaussRF(x=p, y=p, grid=TRUE, model=model, param=param) > get(getOption("device"))(height=4,width=4); > get(getOption("device"))(height=4,width=4); > get(getOption("device"))(height=4,width=4); > > # another grid, where values are to be simulated > step <- 0.25 # or 0.3 > x <- seq(0, 7, step) > > # standardisation of the output > lim <- range( c(x, p) ) > zlim <- c(-2.6, 2.6) > colour <- rainbow(100) > > ## visualise generated spatial data > dev.set(2) postscript 2 > image(p, p, data, xlim=lim, ylim=lim, zlim=zlim, col=colour) > > #conditional simulation > krige.method <- "O" ## random field assumption corresponding to > ## those of ordinary kriging > cz <- CondSimu(krige.method, x, x, grid=TRUE, + model=model, param=param, + given=expand.grid(p,p),# if data are given on a grid + # then expand the grid first + data=data) .................................................................................... > range(cz) [1] -2.217608 2.933812 > dev.set(3) postscript 3 > image(x, x, cz, col=colour, xlim=lim, ylim=lim, zlim=zlim) > > #conditional simulation with error term > cze <- CondSimu(krige.method, x, x, grid=TRUE, + model=model, param=c(0, 1/2, 0, 1), + err.model="gauss", err.param=c(0, 1/2, 0, 1), + given=expand.grid(p,p), + data=data) .................................................................................... > range(cze) [1] -2.280828 3.184791 > dev.set(4) postscript 4 > image(x, x, cze, col=colour, xlim=lim, ylim=lim, zlim=zlim) > > > > cleanEx(); ..nameEx <- "CovarianceFct" > > ### * CovarianceFct > > flush(stderr()); flush(stdout()) > > ### Name: CovarianceFct > ### Title: Covariance And Variogram Models > ### Aliases: CovarianceFct Variogram PrintModelList GetModelList > ### GetModelNames > ### Keywords: spatial > > ### ** Examples > > PrintModelList() List of models ============== [See also PrintMethodList()] Circ cut intr TBM2 TBM3 spec drct nugg add hyp oth bessel X - - - - X X - - - - cauchy X - - X X - X - - - - cauchytbm X - - + X - X - - - - circular X - - + - - X - X - - cone - - - - - - - - X - - cubic X - - + X - X - - - - dampedcosine X - - + X - X - - - - exponential X - - X X X X - - X - FD X - - - - - X - - - - fractalB - - X - - - - - - - - fractgauss X - - - - - X - - - - gauss X - - + X X X - X - - gencauchy X X X + X - X - - - - gengneiting X - - + X - X - - - - gneiting X - - + X - X - - - - hyperbolic X - - + X - X - - - - lgd1 X - - + - - X - - - - nsst X - - X X - X - - - - nsst2 X - - + X - X - - - - nugget X X X penta X - - + X - X - - - - power X - - + X - X - - - - qexponential X - - + X - X - - - - spherical X - - X X - X - X - - stable X X X + X - X - - - - wave X - - - - X X - - - - whittlematern X X X + X X X - - - - Circ cut intr TBM2 TBM3 spec drct nugg add hyp oth Legend:'-': method not available 'X': method available for at least some parameters '+': parts are evaluated only approximatively 'o': given method is ignored and an alternative one is used > x <- 0:100 > > # the following five model definitions are the same! > ## (1) very traditional form > (cv <- CovarianceFct(x, model="bessel", c(NA,2,1,5,0.5))) [1] 3.000000000 1.986693308 1.947091712 1.882141578 1.793390227 [6] 1.682941970 1.553398477 1.407785329 1.249467004 1.082052923 [11] 0.909297427 0.734996731 0.562885984 0.396539517 0.239277250 [16] 0.094080005 -0.036483840 -0.150318295 -0.245844691 -0.322030469 [21] -0.378401248 -0.415036082 -0.432546397 -0.432039567 -0.415068587 [26] -0.383569710 -0.339790252 -0.286209069 -0.225452371 -0.160207648 [31] -0.093138499 -0.026803033 0.036421627 0.094406474 0.145327456 [36] 0.187710457 0.220463296 0.242894080 0.254715703 0.256036755 [41] 0.247339562 0.229446477 0.203475930 0.170790023 0.132935726 [46] 0.091581886 0.048454329 0.005271367 -0.036318079 -0.074791659 [51] -0.108804222 -0.137230331 -0.159197398 -0.174108570 -0.181654857 [56] -0.181816401 -0.174853166 -0.161285706 -0.141866999 -0.117546625 [61] -0.089428820 -0.058726112 -0.026710351 0.005336992 0.036173410 [66] 0.064641083 0.089708108 0.110503864 0.126347326 0.136767488 [71] 0.141515337 0.140567134 0.134119136 0.122574133 0.106520550 [76] 0.086705045 0.063999827 0.039366020 0.013814571 -0.011633779 [81] -0.035987915 -0.058323702 -0.077817888 -0.093777359 -0.105662742 [86] -0.113105587 -0.115918612 -0.114098779 -0.107823238 -0.097438447 [91] -0.083443027 -0.066465145 -0.047235394 -0.026556308 -0.005269749 [96] 0.015776548 0.035761972 0.053924306 0.069588124 0.082189266 [101] 0.091294525 > > ## (2) traditional form in list notation > model <- list(model="bessel", param=c(NA,2,1,5,0.5)) > cv - CovarianceFct(x, model=model) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > ## (3) nested model definition > cv - CovarianceFct(x, model="bessel", + param=cbind(c(2, 5, 0.5), c(1, 0, 0))) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > #### most general notation in form of lists > ## (4) isotropic notation > model <- list(list(model="bessel", var=2, kappa=0.5, scale=5), + "+", + list(model="nugget", var=1)) > cv - CovarianceFct(x, model=model) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > ## (5) anisotropic notation > model <- list(list(model="bessel", var=2, kappa=0.5, aniso=0.2), + "+", + list(model="nugget", var=1, aniso=1)) > cv - CovarianceFct(as.matrix(x), model=model) [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [38] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [75] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 > > # The model gneitingdiff was defined in RandomFields v1.0. > # This isotropic covariance function is valid for dimensions less > # than or equal to 3 and has two positive parameters. > # It is a class of models with compact support that allows for > # smooth parametrisation of the differentiability up to order 6. > # The former model `gneitingdiff' must now be coded as > gneitingdiff <- function(p){ + list(list(m="gneiting", v=p[2], s=p[6]*p[4]*10*sqrt(2)/47), + "*", + list(m="whittle", k=p[5], v=1.0, s=p[4]), + "+", + list(m="nugget", v=p[3])) + } > # and then > param <- c(NA, runif(5,max=10)) ## as usual, here an example > CovarianceFct(x,gneitingdiff(param)) [1] 6.376326e+00 2.443006e+00 1.906663e+00 1.262948e+00 7.048018e-01 [6] 3.248597e-01 1.190231e-01 3.235830e-02 5.746833e-03 5.126394e-04 [11] 1.167053e-05 3.321299e-09 0.000000e+00 0.000000e+00 0.000000e+00 [16] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [21] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [26] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [31] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [36] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [41] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [46] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [51] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [56] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [61] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [66] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [71] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [76] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [81] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [86] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [91] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [96] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 [101] 0.000000e+00 > ## instead of formerly CovarianceFct(x,"gneitingdiff",param) > > > > cleanEx(); ..nameEx <- "Dev" > > ### * Dev > > flush(stderr()); flush(stdout()) > > ### Name: Dev > ### Title: Choosing the device > ### Aliases: Dev > ### Keywords: device utilities > > ### ** Examples > > ## first an eps-file test.eps is created, then a jpeg-file, > ## finally the figure is plotted on the screen > dir(pattern="test*") character(0) > dev.list <- list(TRUE, 1) > if (interactive()) dev.list <- c(dev.list, "jpeg") > for (dev in dev.list) { + print(dev) + size <- if (dev=="jpeg") 450 else 5 + err <- class(try(Dev(TRUE, dev, ps="test", height=size, width=size))) + if (err!="try-error") { + plot(0, 0, main=paste("dev=", dev)) + # readline("press return") + Dev(FALSE) + } + } [1] TRUE creating test.eps [1] 1 > dir(pattern="test*") [1] "test.eps" > > > > cleanEx(); ..nameEx <- "EmpiricalVariogram" > > ### * EmpiricalVariogram > > flush(stderr()); flush(stdout()) > > ### Name: EmpiricalVariogram > ### Title: Empirical (Semi-)Variogram > ### Aliases: EmpiricalVariogram > ### Keywords: spatial > > ### ** Examples > > ############################################################# > ## this example checks whether a certain simulation method ## > ## works well for a specified covariance model and ## > ## a configuration of points ## > ############################################################# > x <- seq(0, 10, 0.5) > y <- seq(0, 10, 0.5) > gridtriple <- FALSE ## see help("GaussRF") > model <- "whittle" ## whittlematern > bins <- seq(0, 5, 0.001) > repetition <- 5 ## by far too small to get reliable results!! > ## It should be of order 500, > ## but then it will take some time > ## to do the simulations > param <- c(mean=1, variance=10, nugget=5, scale=2, alpha=2) > f <- GaussRF(x=x, y=y, grid=TRUE, gridtriple=gridtriple, + model=model, param=param, method="TBM3", + n=repetition) *****> binned <- EmpiricalVariogram(x=x, y=y, data=f, grid=TRUE, + gridtriple=gridtriple, bin=bins) > truevariogram <- Variogram(binned$c, model, param) > matplot(binned$c, cbind(truevariogram,binned$e), pch=c("*","e")) > ##black curve gives the theoretical values > > > > > cleanEx(); ..nameEx <- "GaussRF" > > ### * GaussRF > > flush(stderr()); flush(stdout()) > > ### Name: GaussRF > ### Title: Gaussian Random Fields > ### Aliases: GaussRF InitGaussRF > ### Keywords: spatial > > ### ** Examples > > ############################################################# > ## ## > ## Examples using the symmetric stable model, also called ## > ## "powered exponential model" ## > ## ## > ############################################################# > PrintModelList() ## the complete list of implemented models List of models ============== [See also PrintMethodList()] Circ cut intr TBM2 TBM3 spec drct nugg add hyp oth bessel X - - - - X X - - - - cauchy X - - X X - X - - - - cauchytbm X - - + X - X - - - - circular X - - + - - X - X - - cone - - - - - - - - X - - cubic X - - + X - X - - - - dampedcosine X - - + X - X - - - - exponential X - - X X X X - - X - FD X - - - - - X - - - - fractalB - - X - - - - - - - - fractgauss X - - - - - X - - - - gauss X - - + X X X - X - - gencauchy X X X + X - X - - - - gengneiting X - - + X - X - - - - gneiting X - - + X - X - - - - hyperbolic X - - + X - X - - - - lgd1 X - - + - - X - - - - nsst X - - X X - X - - - - nsst2 X - - + X - X - - - - nugget X X X penta X - - + X - X - - - - power X - - + X - X - - - - qexponential X - - + X - X - - - - spherical X - - X X - X - X - - stable X X X + X - X - - - - wave X - - - - X X - - - - whittlematern X X X + X X X - - - - Circ cut intr TBM2 TBM3 spec drct nugg add hyp oth Legend:'-': method not available 'X': method available for at least some parameters '+': parts are evaluated only approximatively 'o': given method is ignored and an alternative one is used > model <- "stable" > mean <- 0 > variance <- 4 > nugget <- 1 > scale <- 10 > alpha <- 1 ## see help("CovarianceFct") for additional > ## parameters of the covariance functions > step <- 1 ## nicer, but also time consuming if step <- 0.1 > x <- seq(0, 20, step) > y <- seq(0, 20, step) > f <- GaussRF(x=x, y=y, model=model, grid=TRUE, + param=c(mean, variance, nugget, scale, alpha)) > image(x, y, f) > > ############################################################# > ## ... using gridtriple > step <- 1 ## nicer, but also time consuming if step <- 0.1 > x <- c(0, 20, step) ## note: vectors of three values, not a > y <- c(0, 20, step) ## sequence > f <- GaussRF(grid=TRUE, gridtriple=TRUE, + x=x ,y=y, model=model, + param=c(mean, variance, nugget, scale, alpha)) > image(seq(x[1],x[2],x[3]), seq(y[1],y[2],y[3]), f) > > ############################################################# > ## arbitrary points > x <- runif(100, max=20) > y <- runif(100, max=20) > z <- runif(100, max=20) # 100 points in 3 dimensional space > (f <- GaussRF(grid=FALSE, + x=x, y=y, z=z, model=model, + param=c(mean, variance, nugget, scale, alpha))) [1] -1.765530964 -1.074449547 1.087728574 -0.141320633 -0.797834064 [6] 1.437393503 2.843176016 0.681962189 0.474598599 0.681783508 [11] -2.060081566 -1.837517389 3.352547218 1.421692998 0.944839996 [16] 1.538052340 0.129844713 1.119976083 0.460087896 -0.235957527 [21] -1.432894415 2.714452038 2.501969800 0.275179620 1.277719002 [26] -1.340230717 3.341004266 -0.832345004 -3.266687610 2.700571944 [31] 3.330667464 -0.345819039 -1.492902504 -3.113486280 -0.022925489 [36] -2.588247128 -0.694878503 1.712092176 1.898859253 0.413845835 [41] 0.588525947 -0.223560204 3.414581124 -0.699760358 -2.761209347 [46] -3.242620417 0.751303558 0.646014438 1.104851010 0.239166506 [51] 0.792605716 -0.633470293 -0.541786701 -0.921884435 3.313603537 [56] -1.461122355 -0.167225593 0.962070680 1.603602296 1.427070954 [61] -1.411688407 0.108055950 0.167248243 -2.899429751 1.928913585 [66] 2.748026274 -1.774694367 0.272642364 0.684473441 -2.256908537 [71] 2.022902899 -0.008877589 0.088571603 1.082747087 -2.398251968 [76] -1.543498735 0.654289380 1.248651151 -0.462068048 1.667493654 [81] 1.455111818 -0.482944541 -0.148260811 2.083267743 -0.565519278 [86] 0.996526541 -0.824413113 -1.598146039 -0.301511319 -0.856824426 [91] -0.349527433 -1.714731026 2.629213836 -1.974542364 -0.854623140 [96] -4.065433014 -0.210080911 -0.737056677 1.489141985 -1.512859715 > > ############################################################# > ## usage of a specific method > ## -- the complete list can be obtained by PrintMethodList() > x <- runif(100, max=20) # arbitrary points > y <- runif(100, max=20) > (f <- GaussRF(method="dir", # direct matrix decomposition + x=x, y=y, model=model, grid=FALSE, + param=c(mean, variance, nugget, scale, alpha))) [1] 0.56537740 -1.27663110 -1.26782016 1.43293788 -1.34898557 2.18947052 [7] 0.17215652 -1.44757041 -1.93765049 1.64928390 -1.82143676 0.29490603 [13] 0.68895454 0.68598908 1.26359035 -1.66575714 -0.00065955 -1.04377291 [19] -2.40481105 1.97029517 3.34414115 2.02560889 -1.05860027 -0.25710012 [25] -0.48163982 0.08094776 -0.28553967 -3.13897783 -1.58227040 -3.22364610 [31] -2.52622652 -3.64517833 2.10726820 -2.65081584 0.47878144 -2.19065753 [37] 1.20641237 1.09825275 0.89488171 0.17171917 2.14091063 -1.31843030 [43] 5.21560761 2.02422607 0.68412701 1.59821764 1.57015456 -0.10447846 [49] -4.13606921 -2.01041149 -1.88188825 -0.54739106 0.19638657 -1.08933964 [55] -0.06414943 -1.12889804 6.23602882 1.10041774 -0.29236278 -0.94257103 [61] -2.10397317 2.22851688 0.78291171 -1.30026994 -1.04674902 0.46999679 [67] 0.52637536 1.49125337 2.25002552 0.88118318 -1.73799617 -0.41349685 [73] 1.04444127 -2.66569885 -1.29477072 -2.35740796 -2.33421167 -0.16240186 [79] 1.17761381 0.32523718 -1.41397176 2.33433220 -0.66902707 -0.18985613 [85] 0.94638690 2.33595915 0.46060055 3.31476693 1.22834281 1.55901534 [91] 0.85056199 2.40122504 1.92745577 0.91326225 -2.19310876 -0.20590081 [97] -2.02278099 0.61156859 -1.12590960 -0.68512069 > > ############################################################# > ## simulating several random fields at once > step <- 1 ## nicer, but also time consuming if step <- 0.1 > x <- seq(0, 20, step) # grid > y <- seq(0, 20, step) > f <- GaussRF(n=3, # three simulations at once + x=x, y=y, model=model, grid=TRUE, + param=c(mean, variance, nugget, scale, alpha)) ***> image(x, y, f[,,1]) > image(x, y, f[,,2]) > image(x, y, f[,,3]) > > > ############################################################# > ## ## > ## Examples using the extended definition forms ## > ## ## > ## ## > ############################################################# > > ## note that the output seems plausible but not checked!!!! > > ## tbm may also be used for multiplicate models (if they have > ## *exactly* the same anisotropy parameters) > x <- (0:100)/10 > m <- matrix(c(1,2,3,4),ncol=2)/5 > z <- GaussRF(x=x, y=x, grid=TRUE, + model=list( + list(m="power",v=1,k=2,a=m), + "*", list(m="sph", v=1, a=m) + ), + me="TBM3", reg=0,n=1) > print(c(mean(as.double(z)),var(as.double(z)))) [1] -0.1381594 1.0210547 > image(z,zlim=c(-3,3)) > > ## non-separable space-time model applied for two space dimensions > ## note that tbm method does not work nicely, but at least > ## in some special cases. > x <- y <- (1:32)/2 ## grid definition, but as a sequence > T <- c(1,32,1)*10 ## note necessarily gridtriple definition > aniso <- diag(c(0.5,8,1)) > k <- c(1,phi=1,1,1,psi=1,dim=2) > model <- list(list(m="nsst", v=1, k=k, a=aniso)) > z <- GaussRF(x=x, y=y, T=T, grid=TRUE, model=model) > rl <- function() if (interactive()) readline("Press return") > for (i in 1:dim(z)[3]) { image(z[,,i]); rl();} > for (i in 1:dim(z)[2]) { image(z[,i,]); rl();} > for (i in 1:dim(z)[1]) { image(z[i,,]); rl();} > > ############################################################# > ## ## > ## Brownian motion ## > ## using Stein's method ## > ## ## > ############################################################# > # 2d > step <- 0.3 ## nicer, but also time consuming if step <- 0.1 > x <- seq(0, 10, step) > kappa <- 1 # in [0,2) > z <- GaussRF(x=x, y=x, grid=TRUE, model="fractalB", + param=c(0,1,0,1,kappa)) Attention: Intrinsic embedding is to be applied, so the simulated random field will not be stationary. > image(z,zlim=c(-3,3)) > > # 3d > x <- seq(0, 3, step) > kappa <- 1 # in [0,2) > z <- GaussRF(x=x, y=x, z=x, grid=TRUE, model="fractalB", + param=c(0,1,0,1,kappa)) Attention: Intrinsic embedding is to be applied, so the simulated random field will not be stationary. > rl <- function() if (interactive()) readline("Press return") > for (i in 1:dim(z)[1]) { image(z[i,,]); rl();} > > > ############################################################# > ## This example shows the benefits from stored, ## > ## intermediate results: in case of the circulant ## > ## embedding method, the speed is doubled in the second ## > ## simulation. ## > ############################################################# > DeleteAllRegisters() > RFparameters(Storing=TRUE,PrintLevel=1) > y <- x <- seq(0, 50, 0.2) > (p <- c(runif(3), runif(1)+1)) [1] 0.3839134 0.9484392 0.5549689 1.8277623 > ut <- system.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen", + method="circ", param=p)) > image(x, y, f) > hist(f) > c( mean(as.vector(f)), var(as.vector(f)) ) [1] 0.4189124 1.6323272 > cat("system time (first call)", format(ut,dig=3),"\n") system time (first call) 1.12 0.01 1.14 0.00 0.00 > > # second call with the *same* parameters is much faster: > ut <- system.time(f <- GaussRF(x=x,y=y,grid=TRUE,model="exponen", + method="circ",param=p)) > image(x, y, f) > hist(f) > c( mean(as.vector(f)), var(as.vector(f)) ) [1] 0.2348063 1.3859838 > cat("system time (second call)", format(ut,dig=3),"\n") system time (second call) 0.64 0.00 0.66 0.00 0.00 > > > > > > cleanEx(); ..nameEx <- "GetPracticalRange" > > ### * GetPracticalRange > > flush(stderr()); flush(stdout()) > > ### Name: GetPracticalRange > ### Title: Determination of the practical range > ### Aliases: GetPracticalRange PracticalRange > ### Keywords: spatial > > ### ** Examples > > GetPracticalRange("exponential") [1] 2.995732 > GetPracticalRange("whittlematern",kappa=2) [1] 5.475091 > GetPracticalRange("gengneiting",kappa=c(a=3,b=5)) [1] 0.5140105 > > > > cleanEx(); ..nameEx <- "GetRegisterInfo" > > ### * GetRegisterInfo > > flush(stderr()); flush(stdout()) > > ### Name: GetRegisterInfo > ### Title: Internal information > ### Aliases: GetRegisterInfo > ### Keywords: spatial > > ### ** Examples > > str(GetRegisterInfo(0)) List of 13 $ active : logi TRUE $ grid : logi TRUE $ anisotropy : logi FALSE $ time.comp : logi FALSE $ distrib. : chr "Gauss" $ model :List of 2 ..$ :List of 4 .. ..$ name : chr "exponential" .. ..$ method : chr "circulant embedding" .. ..$ param : num [1:10] 0.948 0.000 0.000 0.000 0.000 ... .. ..$ processed: logi TRUE ..$ :List of 5 .. ..$ op : chr "+" .. ..$ name : chr "nugget" .. ..$ method : chr "nugget" .. ..$ param : num [1:10] 0.555 0.000 0.000 0.000 0.000 ... .. ..$ processed: logi TRUE $ mean : num 0.384 $ TBMline : chr "circulant embedding" $ spatial.dim : int 2 $ timespacedim: int 2 $ spatial.pts : int 63001 $ total.pts : int 63001 $ trend.mode : int 0 > try(GaussRF(1:4, grid=TRUE, model="exp", param=c(1,2,3,4))) [1] -0.4007933 0.9538298 -0.9902084 3.7140697 > str(GetRegisterInfo(0)) List of 13 $ active : logi TRUE $ grid : logi TRUE $ anisotropy : logi FALSE $ time.comp : logi FALSE $ distrib. : chr "Gauss" $ model :List of 2 ..$ :List of 4 .. ..$ name : chr "exponential" .. ..$ method : chr "direct matrix decomposition" .. ..$ param : num [1:10] 2 0 0 0 0 0 0 0 4 0.25 .. ..$ processed: logi TRUE ..$ :List of 5 .. ..$ op : chr "+" .. ..$ name : chr "nugget" .. ..$ method : chr "direct matrix decomposition" .. ..$ param : num [1:10] 3 0 0 0 0 0 0 0 1 1 .. ..$ processed: logi TRUE $ mean : num 1 $ TBMline : chr "circulant embedding" $ spatial.dim : int 1 $ timespacedim: int 1 $ spatial.pts : int 4 $ total.pts : int 4 $ trend.mode : int 0 > > try(GaussRF(runif(4), grid=FALSE, model="exp", param=c(1,2,3,4), me="ci")) [1] 1.8442519 2.2788159 1.9972826 0.4710518 > str(GetRegisterInfo(0, TRUE)) List of 13 $ active : logi TRUE $ grid : logi FALSE $ anisotropy : logi FALSE $ time.comp : logi FALSE $ distrib. : chr "Gauss" $ model :List of 2 ..$ :List of 4 .. ..$ name : chr "exponential" .. ..$ method : chr "circulant embedding" .. ..$ param : num [1:10] 2 0 0 0 0 0 0 0 4 0.25 .. ..$ processed: logi FALSE ..$ :List of 5 .. ..$ op : chr "+" .. ..$ name : chr "nugget" .. ..$ method : chr "nugget" .. ..$ param : num [1:10] 3 0 0 0 0 0 0 0 1 1 .. ..$ processed: logi TRUE $ mean : num 1 $ TBMline : chr "circulant embedding" $ spatial.dim : int 1 $ timespacedim: int 1 $ spatial.pts : int 4 $ total.pts : int 4 $ trend.mode : int 0 > > > > cleanEx(); ..nameEx <- "Kriging" > > ### * Kriging > > flush(stderr()); flush(stdout()) > > ### Name: Kriging > ### Title: Kriging methods > ### Aliases: Kriging > ### Keywords: spatial > > ### ** Examples > > ## creating random variables first > ## here, a grid is chosen, but does not matter > step <- 0.25 > x <- seq(0,7,step) > param <- c(0,1,0,1) > model <- "exponential" > RFparameters(PracticalRange=FALSE) > p <- 1:7 > points <- as.matrix(expand.grid(p,p)) > data <- GaussRF(points, grid=FALSE, model=model, param=param) > > ## visualise generated spatial data > zlim <- c(-2.6,2.6) > colour <- rainbow(100) > image(p, p, xlim=range(x), ylim=range(x), + matrix(data,ncol=length(p)), + col=colour,zlim=zlim) > > ## now: kriging > krige.method <- "O" ## ordinary kriging > z <- Kriging(krige.method=krige.method, + x=x, y=x, grid=TRUE, + model=model, param=param, + given=points, data=data) .................................................................................... > image(x,x,z,col=colour,zlim=zlim) > > > > cleanEx(); ..nameEx <- "Locator" > > ### * Locator > > flush(stderr()); flush(stdout()) > > ### Name: Locator > ### Title: Graphical Input > ### Aliases: Locator > ### Keywords: iplot utilities > > ### ** Examples > > ## see useraction > > > > cleanEx(); ..nameEx <- "MaxStableRF" > > ### * MaxStableRF > > flush(stderr()); flush(stdout()) > > ### Name: MaxStableRF > ### Title: Max-Stable Random Fields > ### Aliases: MaxStableRF InitMaxStableRF > ### Keywords: spatial > > ### ** Examples > > n <- 30 ## nicer, but time consuming if n <- 100 > x <- y <- 1:n > ms0 <- MaxStableRF(x, y, grid=TRUE, model="exponen", + param=c(0,1,0,40), maxstable="extr" + ,CE.force = TRUE ) > image(x,y,ms0) > > ############################################################ > ## ## > ## Plots used in M. Schlather, Extremes 5:1, 33-44, 2002 ## > ## ## > ############################################################ > > pts <- if (interactive()) 512 else 32 > x <- (1:pts) / pts * 10 > scalegauss <- 1.5 > RFparameters(MPP.radius=2*scalegauss, Print=3) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > runif(1) [1] 0.2873062 > save.seed <- .Random.seed > ms1 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, scalegauss), + maxstable="Bool", grid=TRUE) Note: area has been enlarged by fixed value (3.000000)return code 0 : Message: fine. 100: 0 0.002585 0.009122 200: 0 0.004149 0.005213 300: 7 0.003379 0.003461 400: 32 0.002403 0.002627 WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > image(x, x, sqrt(ms1)) > > scalecone <- 3.7 > .Random.seed <- save.seed > ms2 <- MaxStableRF(x, x, model="cone", param=c(0, 1, 0, scalecone, + 0, 0, 1), maxstable="Bool", grid=TRUE) Note: area has been enlarged by fixed value (3.000000)return code 0 : Message: fine. 100: 0 0.003257 0.009122 200: 0 0.004257 0.005213 300: 32 0.002954 0.003461 WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > image(x, x, sqrt(ms2)) > > .Random.seed <- save.seed > ms3 <- MaxStableRF(x, x, model="exponen", param=c(0, 1, 0, 1), + maxstable="extr", grid=TRUE, method="ci") m[0]=64, m[1]=64, mtot=4096 return code 0 : Message: fine. 10: 22 8.919853e-02 2.675956e-01 20: 22 5.144433e-02 1.543330e-01 30: 263 3.092368e-02 9.277103e-02 40: 836 2.245762e-02 6.737287e-02 WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > image(x, x, sqrt(ms3)) > > .Random.seed <- save.seed > ms4 <- MaxStableRF(x, x, model="gauss", param=c(0, 1, 0, 1), + maxstable="extr", grid=TRUE, method="ci") m[0]=64, m[1]=64, mtot=4096 return code 0 : Message: fine. 10: 19 8.919853e-02 2.675956e-01 20: 677 5.144433e-02 1.543330e-01 30: 772 3.092368e-02 9.277103e-02 WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > image(x, x, sqrt(ms4)) > > > > > cleanEx(); ..nameEx <- "RFMethods" > > ### * RFMethods > > flush(stderr()); flush(stdout()) > > ### Name: RFMethods > ### Title: Simulation Techniques > ### Aliases: RFMethods PrintMethodList GetMethodNames > ### Keywords: spatial > > ### ** Examples > > PrintMethodList() Methods for generating Gaussian random fields ============================================= * circulant embedding * cutoff CE * intrinsic CE * TBM2 (2dim. turning bands) * TBM3 * spectral TBM (2dim) * direct matrix decomposition * nugget * add.MPP * hyperplanes * other methods(not implem.) Methods for non-Gaussian random fields ====================================== * max.MPP > > > > cleanEx(); ..nameEx <- "RFparameters" > > ### * RFparameters > > flush(stderr()); flush(stdout()) > > ### Name: RFparameters > ### Title: Control Parameters > ### Aliases: RFparameters > ### Keywords: spatial > > ### ** Examples > > RFparameters(Storing=TRUE) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > str(RFparameters()) List of 61 $ PracticalRange : logi FALSE $ PrintLevel : int 3 $ pch : chr "*" $ Storing : logi TRUE $ stationary.only : logi NA $ exactness : logi NA $ CE.force : logi FALSE $ CE.mmin : int [1:4] 0 0 0 0 $ CE.strategy : int 0 $ CE.maxmem : num 2e+07 $ CE.tolIm : num 0.001 $ CE.tolRe : num -1e-07 $ CE.trials : int 3 $ CE.userfft : logi FALSE $ cutoff.a : num 0 $ intrinsic.r : num 0 $ direct.checkprecision : logi FALSE $ direct.bestvariables : int 600 $ direct.maxvariables : int 1800 $ direct.method : int 0 $ direct.requiredprecision: num 1e-11 $ spectral.grid : logi TRUE $ spectral.lines : int 500 $ TBM.method : chr "circulant embedding" $ TBM2.every : int 0 $ TBM2.lines : int 60 $ TBM2.linesimufactor : num 2 $ TBM2.linesimustep : num 0 $ TBM2.num : logi TRUE $ TBM3D2.every : int 0 $ TBM3D2.lines : int 500 $ TBM3D2.linesimufactor : num 2 $ TBM3D2.linesimustep : num 0 $ TBM3D3.every : int 0 $ TBM3D3.lines : int 500 $ TBM3D3.linesimufactor : num 2 $ TBM3D3.linesimustep : num 0 $ TBMCE.force : logi FALSE $ TBMCE.mmin : int [1:4] 0 0 0 0 $ TBMCE.strategy : int 0 $ TBMCE.maxmem : num 1e+07 $ TBMCE.tolIm : num 0.001 $ TBMCE.tolRe : num -1e-07 $ TBMCE.trials : int 3 $ TBMCE.userfft : logi TRUE $ add.MPP.realisations : num 100 $ MPP.approxzero : num 0.001 $ MPP.radius : num 3 $ hyper.superpos : int 300 $ hyper.maxlines : int 1000 $ hyper.mar.distr : int 0 $ hyper.mar.param : num 0 $ maxstable.maxGauss : num 3 $ covmaxchar : int 14 $ covnr : int 27 $ distrmaxchar : int 10 $ distrnr : int 3 $ maxdim : int 4 $ maxmodels : int 10 $ methodmaxchar : int 31 $ methodnr : int 11 > > > > cleanEx(); ..nameEx <- "Readline" > > ### * Readline > > flush(stderr()); flush(stdout()) > > ### Name: Readline > ### Title: Read a Line > ### Aliases: Readline > ### Keywords: utilities > > ### ** Examples > > ## see useraction > > > > cleanEx(); ..nameEx <- "Registers" > > ### * Registers > > flush(stderr()); flush(stdout()) > > ### Name: DeleteRegister > ### Title: Deleting Intermediate Results > ### Aliases: DeleteRegister DeleteAllRegisters > ### Keywords: spatial > > ### ** Examples > > DeleteAllRegisters() ## no effect visible for the user > > > > cleanEx(); ..nameEx <- "ShowModels" > > ### * ShowModels > > flush(stderr()); flush(stdout()) > > ### Name: ShowModels > ### Title: Interactive Choice of Models and Parameters > ### Aliases: ShowModels > ### Keywords: spatial > > ### ** Examples > > # first example: one-dimensional simulations > > options(locatorBell=FALSE) > x <- seq(1,10,0.1); > ShowModels(x=x) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. NULL > > # second example: two-dimensional simulations and > # empirical variogram > dx <- runif(300,0,8) > dy <- runif(300,0,8) > dz <- GaussRF(x=dx, y=dy, grid=FALSE, model="gaus", + param=c(1,2,1,2)) method=Cholesky return code 0 : Message: fine. WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > ev <- EmpiricalVariogram(x=dx, y=dy, data=dz, grid=FALSE, + bin=(-1:20)/4) > x <- seq(1,5,0.1); > ShowModels(x=x, y=x, empirical=ev) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. NULL > > # third example: two-dimensional anistropic simulations and > # link function > x <- seq(1,10,0.1) > ShowModels(x=x, y=x, link=function(x) exp(x), + model=list(list(model="spheric", var=1, aniso=c(1,0,0,5)))) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. $model $model[[1]] $model[[1]]$model [1] "spherical" $model[[1]]$var [1] 1 $model[[1]]$kappas NULL $model[[1]]$aniso [1] 1 0 0 5 $model$op [1] "+" $model[[3]] $model[[3]]$model [1] "nugget" $model[[3]]$var [1] 0 $model[[3]]$kappas NULL $model[[3]]$aniso [1] 1 0 0 1 $mean [1] 0 $method NULL $trend NULL > > x <- seq(1,10,0.1) > ShowModels(x=x, link=function(x) exp(x), + model=list(list(model="spheric",var=1, scale=1))) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. $model [1] "spherical" $param [1] 0 1 0 1 $method NULL $trend NULL > > x <- seq(1,10,0.1) > ShowModels(x=x, link="MaxStable", fixed.rs=TRUE, + model=list(list(model="gauss",var=1, scale=1)), type="l") Error code F77_CALL(chol) = 9 WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. $model [1] "gauss" $param [1] 0 1 0 1 $method NULL $trend NULL > > > > cleanEx(); ..nameEx <- "auxiliary.fcts" > > ### * auxiliary.fcts > > flush(stderr()); flush(stdout()) > > ### Name: CheckXT > ### Title: Internal functions - do not use them directly > ### Aliases: CheckXT PrepareModel convert.to.readable plotWithCircles > ### GetDistributionNames > ### Keywords: spatial internal > > ### ** Examples > > > x <- function(...) { + str(PrepareModel(...)) + cat("--------------------------------\n") + str(convert.to.readable(PrepareModel(...))) + } > model <- list(list(model="whi", kappa=5, var=2, s=4), "+", + list(model="whi", kappa=1, var=3, s=0)) ## s=0 should not be used only in > ## a model definition where the parameters are > ## are given in a matrix, see the result > x(model=model, ti=1, me="ci") List of 8 $ covnr : int [1:2] 26 26 $ param : num [1:6] 2 5 4 3 1 0 $ anisotropy : logi FALSE $ op : num 0 $ mean : num 0 $ trend : NULL $ method : int [1:2] 0 0 $ timespacedim: int 1 -------------------------------- List of 5 $ model : chr "whittlematern" $ param : num [1:3, 1:2] 2 4 5 3 0 1 $ method: chr "circulant embedding" $ mean : num 0 $ trend : NULL > > ## since convert.to.readable performs a one-step simplification, > ## iterative calls may further simplify the model > xx <- convert.to.readable(PrepareModel(model=model, ti=1, me="ci")) > x(model=xx$mo, pa=xx$pa, ti=1, me=xx$me) List of 8 $ covnr : int [1:2] 26 19 $ param : num [1:5] 2 5 4 3 1 $ anisotropy : logi FALSE $ op : num 0 $ mean : num 0 $ trend : NULL $ method : num -1 $ timespacedim: int 1 -------------------------------- List of 4 $ model : chr "whittlematern" $ param : num [1:5] 0 2 3 4 5 $ method: NULL $ trend : NULL > > ## back to the matrix definition of nested models > str(convert.to.readable(PrepareModel(xx, ti=1), allowed="nested")) nugget simulation method overwritten. List of 5 $ model : chr "whittlematern" $ param : num [1:3, 1:2] 2 4 5 3 0 0 $ method: chr "circulant embedding" $ mean : num 0 $ trend : NULL > > ## back to the (correct) list definition > str(convert.to.readable(PrepareModel(xx, ti=1), allowed="list")) List of 4 $ model :List of 3 ..$ :List of 4 .. ..$ model : chr "whittlematern" .. ..$ var : num 2 .. ..$ kappas: num 5 .. ..$ scale : num 4 ..$ op: chr "+" ..$ :List of 4 .. ..$ model : chr "nugget" .. ..$ var : num 3 .. ..$ kappas: NULL .. ..$ scale : num 1 $ mean : num 0 $ method: chr [1:2] "circulant embedding" "circulant embedding" $ trend : NULL > > > > cleanEx(); ..nameEx <- "eval.parameters" > > ### * eval.parameters > > flush(stderr()); flush(stdout()) > > ### Name: eval.parameters > ### Title: Interactive menu > ### Aliases: eval.parameters > ### Keywords: spatial > > ### ** Examples > > ## the following lines define a menu that does not make > ## too much sense, but shows the various kinds of buttons > > quadratic <- function(d, v, a, mini=0, maxi=Inf) { + d <- pmin(1, pmax(0, d)) - 0.5 + d <- ((d>0) * 2 - 1) * d^2 * a * 4 + if (missing(v)) d else pmax(mini, pmin(maxi, v + d)) + } > > simulate <- function(H, par) { ## not a serious example + print(c(H$x$var, par, runif(1))) + return(H) ## the function must return the first parameter + } > > entry <- list( + list(name="Nonsense Menu"), + list(name="Simulate!", val="simulate", col="blue"), + list(name="show H", val="str(H)", col="blue"), + list(name="colx", var="colour", + val=c("red", "green", "blue", "brown")), + list(name="open", var="closed", val=FALSE, par=4.5), + list(name="modifying", var="modify", val=TRUE, par=5), + list(name="probability", var="probab", delta=FALSE, + val=function(d, v) pmin(1, pmax(0, d))), + list(name="variance", var="var", delta=TRUE, + val=function(d, v) quadratic(d, v, 10)), + list(name="name", var="name", par=3, cond="modify"), + ) > > scr <- split.screen(rbind(c(0, 0.45, 0, 1), c(0.5, 1, 0, 1))) > ## before proceeding make sure that both the screen and the xterm > ## are completely visible > > H <- list(modify=5, x=list()) # note that in this example eval.parameters > ## will be called by by H$x, hence modify=5 will be left > ## unchanged. > options(locatorBell=FALSE) > useraction("start.register") ## registring the user's input > str(eval.parameters("H$x", entry, simulate, update=TRUE, dev=scr[2], + H=H, par=17)) # do not forget to call by name Creating `H$x$name' ... Creating `H$x$var' ... Creating `H$x$probab' ... Creating `H$x$modify' ... Creating `H$x$closed' ... Creating `H$x$colour' ... List of 2 $ modify: num 5 $ x :List of 7 ..$ name : chr "" ..$ var : num 0 ..$ probab : num 0.5 ..$ modify : logi TRUE ..$ closed : logi TRUE ..$ colour : num 1 ..$ .history: list() > getactions() $mismatch [1] "user input replay mismatch" $sleep [1] 2.5 $wait [1] 0.1 $mode [1] "continue.register" $actions $actions[[1]] $actions[[1]]$info [1] "H$x" $n [1] Inf $print [1] 0 > > ## replay the user's input > useraction("replay") > str(eval.parameters("H$x", entry, simulate, update=TRUE, dev=scr[2], + H=H, par=17)) Creating `H$x$name' ... Creating `H$x$var' ... Creating `H$x$probab' ... Creating `H$x$modify' ... Creating `H$x$closed' ... Creating `H$x$colour' ... List of 2 $ modify: num 5 $ x :List of 7 ..$ name : chr "" ..$ var : num 0 ..$ probab : num 0.5 ..$ modify : logi TRUE ..$ closed : logi TRUE ..$ colour : num 1 ..$ .history: list() > > > > cleanEx(); ..nameEx <- "fitvario" > > ### * fitvario > > flush(stderr()); flush(stdout()) > > ### Name: fitvario > ### Title: LSQ and Maximum Likelihood Estimation of Random Field Parameters > ### Aliases: fitvario mleRF fitvario.default > ### Keywords: spatial > > ### ** Examples > > model <-"gencauchy" > param <- c(0, 1, 0, 1, 1, 2) > estparam <- c(0, NA, 0, NA, NA, 2) ## NA means: "to be estimated" > ## sequence in `estparam' is > ## mean, variance, nugget, scale, (+ further model parameters) > ## So, mean, variance, and scale will be estimated here. > ## Nugget is fixed and equals zero. > points <- 100 > x <- runif(points,0,3) > y <- runif(points,0,3) ## 100 random points in square [0, 3]^2 > d <- GaussRF(x=x, y=y, grid=FALSE, model=model, param=param, n=10) method=Cholesky return code 0 : Message: fine. ********** WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, + lower=c(0.1, 0.1), upper=c(1.9, 5), cross.me=NULL)) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. self plain sqrt.nr sd.inv internal1 internal2 internal3 ml reml WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. List of 3 $ ev :List of 7 ..$ centers : num [1:21] 0.0000 0.0459 0.1377 0.2295 0.3212 ... ..$ emp.vario: num [1:21] 0.000 0.120 0.263 0.388 0.503 ... ..$ sd : num [1:21] 0.000 0.171 0.359 0.554 0.671 ... ..$ n.bin : int [1:21] 1000 280 780 1540 2140 2720 2420 3180 3440 3960 ... ..$ phi : NULL ..$ theta : NULL ..$ T : NULL $ variogram:List of 10 ..$ autostart:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.000 0.000 0.217 1.000 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ self :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.174 0.000 1.442 0.883 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ plain :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.00 1.12 0.00 1.07 0.96 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sqrt.nr :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.101 0.000 0.966 0.996 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sd.inv :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.134 0.000 1.226 0.917 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal1:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.151 0.000 1.230 0.928 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal2:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.194 0.000 1.482 0.886 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal3:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.160 0.000 1.301 0.912 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ ml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.934 0.000 0.590 1.130 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ reml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.934 0.000 0.590 1.130 ... .. ..$ method: NULL .. ..$ trend : NULL $ values :List of 10 ..$ autostart: num(0) ..$ self : num 0.0268 ..$ plain : num 0.0173 ..$ sqrt.nr : num 1.12 ..$ sd.inv : num 0.0136 ..$ internal1: num 2.59 ..$ internal2: num 0.105 ..$ internal3: num 0.0403 ..$ ml : num -1020 ..$ reml : num -1020 > > ## The next two estimations give about the same result. > ## For the first the sill is fixed to 1.5. For the second the sill > ## is reached if the estimated variance is smaller than 1.5 > estparam <- c(0, NA, NA, NA, NA, NA) > str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, + sill=1.5, cross.me=NULL)) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. self plain sqrt.nr sd.inv internal1 internal2 internal3 mlError in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 25 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2379758015127500 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 25 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 25 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2379758015127500 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 25 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 25 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2379758015127500 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 25 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 19 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2389758015127503 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 19 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 23 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2369758015127497 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 23 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 25 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2379758015127500 [4] 36.7131962926096804 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 25 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 24 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 1.5000000000000000 2.0000000000000000 7.2379758015127500 [4] 36.7127118081286028 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 24 is not positive de"| __truncated__ ++ reml WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. List of 3 $ ev :List of 7 ..$ centers : num [1:21] 0.0000 0.0459 0.1377 0.2295 0.3212 ... ..$ emp.vario: num [1:21] 0.000 0.120 0.263 0.388 0.503 ... ..$ sd : num [1:21] 0.000 0.171 0.359 0.554 0.671 ... ..$ n.bin : int [1:21] 1000 280 780 1540 2140 2720 2420 3180 3440 3960 ... ..$ phi : NULL ..$ theta : NULL ..$ T : NULL $ variogram:List of 10 ..$ autostart:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.750 0.750 0.573 1.025 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ self :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.4866 0.0134 0.4731 0.9544 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ plain :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.537 0.885 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sqrt.nr :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.562 0.864 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sd.inv :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.520 0.902 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal1:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.503 0.904 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal2:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.00000 1.49288 0.00712 0.45970 0.95221 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal3:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.483 0.923 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ ml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.385 1.103 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ reml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 0.385 1.103 ... .. ..$ method: NULL .. ..$ trend : NULL $ values :List of 10 ..$ autostart: num(0) ..$ self : num 0.0298 ..$ plain : num 0.0213 ..$ sqrt.nr : num 1.39 ..$ sd.inv : num 0.0167 ..$ internal1: num 2.84 ..$ internal2: num 0.108 ..$ internal3: num 0.044 ..$ ml : num -1026 ..$ reml : num -1026 > > estparam <- c(0, NA, NaN, NA, NA, NA) > parampositions(model=model, param=estparam) model: gencauchy param: NA 1 5 4 2 3 > f <- function(param) { + param[5] <- max(0, 1.5 - param[1]) + return(param) + } > str(fitvario(x=cbind(x,y), data=d, model=model, param=estparam, + sill=1, transform=f, cross.me=NULL)) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. self plain sqrt.nr sd.inv internal1LSQ WARNING! forbidden values ! 9.41991 0.8579371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.420381 0.8579371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.419438 0.8579371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.41991 0.8589371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.41991 0.8569371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.41991 0.8579371 0.05 36.71320 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 9.41991 0.8579371 0.05 36.71271 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] internal2 internal3LSQ WARNING! forbidden values ! 0.00943462 0.5680317 0.6389213 3.086472 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5690317 0.6389213 3.086472 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5670317 0.6389213 3.086472 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5680317 0.6399213 3.086472 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5680317 0.6379213 3.086472 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5680317 0.6389213 3.086956 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] LSQ WARNING! forbidden values ! 0.00943462 0.5680317 0.6389213 3.085987 [ 0.00943462 0.05 0.05 0.005844463 , 9.43462 2 10 36.71320 ] mlError in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 18 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.223618610608700 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 18 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 18 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.223618610608700 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 18 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 17 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434149082991759 2.000000000000000 0.050000000000000 31.223618610608700 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 17 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 18 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.223618610608700 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 18 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 17 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032461666 2.000000000000000000 0.051000000000000004 [4] 31.223618610608699697 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 17 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 18 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.223618610608700 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 18 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 17 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.224103095089774 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 17 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 17 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 9.434620814032462 2.000000000000000 0.050000000000000 31.223134126127622 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 17 is not positive de"| __truncated__ reml WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. List of 3 $ ev :List of 7 ..$ centers : num [1:21] 0.0000 0.0459 0.1377 0.2295 0.3212 ... ..$ emp.vario: num [1:21] 0.000 0.120 0.263 0.388 0.503 ... ..$ sd : num [1:21] 0.000 0.171 0.359 0.554 0.671 ... ..$ n.bin : int [1:21] 1000 280 780 1540 2140 2720 2420 3180 3440 3960 ... ..$ phi : NULL ..$ theta : NULL ..$ T : NULL $ variogram:List of 10 ..$ autostart:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.472 1.028 0.484 1.025 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ self :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 2.819 0.775 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ plain :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 3.344 0.699 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sqrt.nr :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 5.293 0.656 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ sd.inv :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.4285 0.0715 0.1876 1.4833 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal1:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 3.044 0.000 32.469 0.527 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal2:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 4.687 0.728 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ internal3:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.500 0.000 3.866 0.718 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ ml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.4285 0.0715 0.1876 1.4833 ... .. ..$ method: NULL .. ..$ trend : NULL ..$ reml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.4285 0.0715 0.1876 1.4833 ... .. ..$ method: NULL .. ..$ trend : NULL $ values :List of 10 ..$ autostart: num(0) ..$ self : num 0.0497 ..$ plain : num 0.0273 ..$ sqrt.nr : num 1.69 ..$ sd.inv : num 0.0139 ..$ internal1: num 5.34 ..$ internal2: num 0.150 ..$ internal3: num 0.0595 ..$ ml : num -1027 ..$ reml : num -1027 > > ## the next call gives a warning, since the user may programme > ## strange things in this setup, and the program cannot check it. > estparam <- c(0, NA, NA, NA, NA, NaN) > parampositions(model=model, param=estparam) model: gencauchy param: NA 1 5 4 2 3 > f <- function(param) {param[3] <- param[2]; param} > system.time(str(fitvario(x=cbind(x,y), data=d, model=model, + param=estparam, transform=f, standard.style=TRUE, + cross.me=NULL), vec.len=6)) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. Warning in fitvario.default(x = x, y = y, z = z, T = T, data = data, model = model, : standard.style and transform!=NULL may cause strange effects -- internal transformation is performed before the user's one self plain sqrt.nr sd.inv internal1 internal2 internal3 mlError in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 47 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 2.000000000000000 34.001556995365405 0.000000000000000 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 47 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 47 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 2.000000000000000 34.001556995365405 0.000000000000000 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 47 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 46 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 2.000000000000000 34.002041479846483 0.000000000000000 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 46 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 40 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 2.00000000000000 34.00107251088433 0.00000000000000 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 40 is not positive de"| __truncated__ Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 47 is not positive definite error in cholesky decomposition -- matrix pos def?[1] 2.000000000000000 34.001556995365405 0.000000000000000 Class 'try-error' chr "Error in chol(CC <- matrix(.C(\"CovarianceMatrixNatSc\", distances, lc, : \n\tthe leading minor of order 47 is not positive de"| __truncated__ reml WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. List of 3 $ ev :List of 7 ..$ centers : num [1:21] 0.0000 0.0459 0.1377 0.2295 0.3212 0.4130 0.5048 0.5966 ... ..$ emp.vario: num [1:21] 0.000 0.120 0.263 0.388 0.503 0.571 0.609 0.643 ... ..$ sd : num [1:21] 0.000 0.171 0.359 0.554 0.671 0.780 0.803 0.896 ... ..$ n.bin : int [1:21] 1000 280 780 1540 2140 2720 2420 3180 3440 3960 3760 4500 4300 4700 4420 ... ..$ phi : NULL ..$ theta : NULL ..$ T : NULL $ variogram:List of 10 ..$ autostart:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 0.5000 0.5000 0.0274 1.0250 1.0250 .. ..$ method: NULL .. ..$ trend : NULL ..$ self :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.126 0.048 0.478 1.148 1.148 .. ..$ method: NULL .. ..$ trend : NULL ..$ plain :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.0783 0.0639 0.4673 1.2193 1.2193 .. ..$ method: NULL .. ..$ trend : NULL ..$ sqrt.nr :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.0326 0.0876 0.4700 1.2990 1.2990 .. ..$ method: NULL .. ..$ trend : NULL ..$ sd.inv :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.1078 0.0503 0.4695 1.1618 1.1618 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal1:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.1540 0.0523 0.5116 1.1236 1.1236 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal2:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.2641 0.0255 0.5632 1.0113 1.0113 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal3:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.1784 0.0415 0.5140 1.0919 1.0919 .. ..$ method: NULL .. ..$ trend : NULL ..$ ml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.0455 0.0475 0.4695 1.1618 1.1618 .. ..$ method: NULL .. ..$ trend : NULL ..$ reml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.0000 1.0455 0.0475 0.4695 1.1618 1.1618 .. ..$ method: NULL .. ..$ trend : NULL $ values :List of 10 ..$ autostart: num(0) ..$ self : num 0.0243 ..$ plain : num 0.0173 ..$ sqrt.nr : num 1.13 ..$ sd.inv : num 0.0132 ..$ internal1: num 2.61 ..$ internal2: num 0.106 ..$ internal3: num 0.0404 ..$ ml : num -1024 ..$ reml : num -1024 [1] 1.23 0.08 1.31 0.00 0.00 > > ## much better programmed, but also much slower: > estmodel <- list(list(model="gencauchy", var=NA, scale=NA, + kappa=list(NA, function(m) m[[1]]$kappa[1]))) > system.time(str(fitvario(x=cbind(x,y), data=d, model=estmodel, + cross.me=NULL), vec.len=6)) WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. self plain sqrt.nr sd.inv internal1 internal2 internal3 ml reml WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. List of 3 $ ev :List of 7 ..$ centers : num [1:21] 0.0000 0.0459 0.1377 0.2295 0.3212 0.4130 0.5048 0.5966 ... ..$ emp.vario: num [1:21] 0.000 0.120 0.263 0.388 0.503 0.571 0.609 0.643 ... ..$ sd : num [1:21] 0.000 0.171 0.359 0.554 0.671 0.780 0.803 0.896 ... ..$ n.bin : int [1:21] 1000 280 780 1540 2140 2720 2420 3180 3440 3960 3760 4500 4300 4700 4420 ... ..$ phi : NULL ..$ theta : NULL ..$ T : NULL $ variogram:List of 10 ..$ autostart:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.943 0.000 0.484 1.025 1.025 .. ..$ method: NULL .. ..$ trend : NULL ..$ self :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.350 0.000 0.601 0.915 0.915 .. ..$ method: NULL .. ..$ trend : NULL ..$ plain :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.218 0.000 0.474 1.030 1.030 .. ..$ method: NULL .. ..$ trend : NULL ..$ sqrt.nr :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.186 0.000 0.449 1.084 1.084 .. ..$ method: NULL .. ..$ trend : NULL ..$ sd.inv :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.269 0.000 0.521 0.965 0.965 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal1:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.269 0.000 0.519 0.997 0.997 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal2:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.353 0.000 0.600 0.934 0.934 .. ..$ method: NULL .. ..$ trend : NULL ..$ internal3:List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 1.295 0.000 0.543 0.969 0.969 .. ..$ method: NULL .. ..$ trend : NULL ..$ ml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.970 0.000 0.326 1.216 1.216 .. ..$ method: NULL .. ..$ trend : NULL ..$ reml :List of 4 .. ..$ model : chr "gencauchy" .. ..$ param : num [1:6] 0.000 0.970 0.000 0.326 1.216 1.216 .. ..$ method: NULL .. ..$ trend : NULL $ values :List of 10 ..$ autostart: num(0) ..$ self : num 0.0289 ..$ plain : num 0.0181 ..$ sqrt.nr : num 1.17 ..$ sd.inv : num 0.0145 ..$ internal1: num 2.66 ..$ internal2: num 0.107 ..$ internal3: num 0.0415 ..$ ml : num -1020 ..$ reml : num -1020 [1] 15.04 0.30 15.47 0.00 0.00 > > > > > cleanEx(); ..nameEx <- "getactions" > > ### * getactions > > flush(stderr()); flush(stdout()) > > ### Name: getactions > ### Title: Get input behaviour > ### Aliases: getactions getactionlist > ### Keywords: utilities > > ### ** Examples > > ## see useraction > > > > cleanEx(); ..nameEx <- "parameter.range" > > ### * parameter.range > > flush(stderr()); flush(stdout()) > > ### Name: parameter.range > ### Title: Range of model specific parameters > ### Aliases: parameter.range > ### Keywords: spatial > > ### ** Examples > > parameter.range("exponential", 1) # NULL NULL > str(parameter.range("power", 1)) List of 2 $ theoretical:List of 1 ..$ : num [1:2, 1] 1 Inf $ practical :List of 1 ..$ : num [1:2, 1] 1 20 > str(parameter.range("power", 2)) List of 2 $ theoretical:List of 1 ..$ : num [1:2, 1] 1.5 Inf $ practical :List of 1 ..$ : num [1:2, 1] 1.5 20 > str(parameter.range("gengneiting", 2)) List of 2 $ theoretical:List of 3 ..$ : num [1:2, 1:2] 1.0 1.0 2.5 Inf ..$ : num [1:2, 1:2] 2.0 2.0 3.5 Inf ..$ : num [1:2, 1:2] 3.0 3.0 4.5 Inf $ practical :List of 3 ..$ : num [1:2, 1:2] 1 1 2.5 20 ..$ : num [1:2, 1:2] 2 2 3.5 20 ..$ : num [1:2, 1:2] 3 3 4.5 20 > str(parameter.range("nsst",1)) List of 2 $ theoretical:List of 9 ..$ : num [1:2, 1:6] 0 2 1 1 0 2 0 1 1 1 ... ..$ : num [1:2, 1:6] 0 Inf 2 2 0 ... ..$ : num [1:2, 1:6] 0 Inf 3 3 0 ... ..$ : num [1:2, 1:6] 0 2 1 1 0 2 0 1 2 2 ... ..$ : num [1:2, 1:6] 0 Inf 2 2 0 ... ..$ : num [1:2, 1:6] 0 Inf 3 3 0 ... ..$ : num [1:2, 1:6] 0 2 1 1 0 2 0 1 3 3 ... ..$ : num [1:2, 1:6] 0 Inf 2 2 0 ... ..$ : num [1:2, 1:6] 0 Inf 3 3 0 ... $ practical :List of 9 ..$ : num [1:2, 1:6] 6e-02 2e+00 1e+00 1e+00 1e-20 ... ..$ : num [1:2, 1:6] 1e-02 1e+01 2e+00 2e+00 1e-20 ... ..$ : num [1:2, 1:6] 9e-02 1e+01 3e+00 3e+00 1e-20 ... ..$ : num [1:2, 1:6] 6e-02 2e+00 1e+00 1e+00 1e-20 ... ..$ : num [1:2, 1:6] 1e-02 1e+01 2e+00 2e+00 1e-20 ... ..$ : num [1:2, 1:6] 9e-02 1e+01 3e+00 3e+00 1e-20 ... ..$ : num [1:2, 1:6] 6e-02 2e+00 1e+00 1e+00 1e-20 ... ..$ : num [1:2, 1:6] 1e-02 1e+01 2e+00 2e+00 1e-20 ... ..$ : num [1:2, 1:6] 9e-02 1e+01 3e+00 3e+00 1e-20 ... > str(parameter.range("nsst2",1)) List of 2 $ theoretical:List of 3 ..$ : num [1:2, 1:7] 0 2 0 Inf 1 ... ..$ : num [1:2, 1:7] 0 2 0 Inf 1 ... ..$ : num [1:2, 1:7] 0 2 0 Inf 1 ... $ practical :List of 3 ..$ : num [1:2, 1:7] 0.05 2.00 0.05 10.00 1.00 ... ..$ : num [1:2, 1:7] 0.05 2.00 0.05 10.00 1.00 ... ..$ : num [1:2, 1:7] 0.05 2.00 0.05 10.00 1.00 ... > str(parameter.range("hyper",1)) List of 2 $ theoretical:List of 1 ..$ : num [1:2, 1:3] 0 Inf -Inf Inf 0 ... $ practical :List of 1 ..$ : num [1:2, 1:3] 1e-06 10 -20 20 1e-06 10 > > > > cleanEx(); ..nameEx <- "parampositions" > > ### * parampositions > > flush(stderr()); flush(stdout()) > > ### Name: parampositions > ### Title: Position of the parameters > ### Aliases: parampositions > ### Keywords: spatial > > ### ** Examples > > parampositions(model="exp", param=c(0,1,0,NA)) model: exponential param: NA 1 NA 2 > ## that is, the nugget in the standard model is removed if naught! > ## the values of the other parameters do not matter. > > parampositions(model="exp", param=c(0,1,NA,NA)) model: exponential param: NA 1 3 2 > parampositions(model="exp", param=c(0,0,1,NA)) model: exponential param: NA 1 3 2 > > parampositions(model="whi", param=cbind(c(1, 1, 1), c(2, 2, 2))) model: whittlematern param: [,1] [,2] [1,] 1 4 [2,] 3 6 [3,] 2 5 > parampositions(model="whi", param=cbind(c(1, 1, 1), c(2, 0, 2))) model: whittlematern param: [,1] [,2] [1,] 1 4 [2,] 3 NA [3,] 2 NA > ## second lines, second column defines a nugget effect since scale is 0! > > try(parampositions(model="whi", param=cbind(c(1, 1, 1), c(0, 0, 0)))) Error in parampositions(model = "whi", param = cbind(c(1, 1, 1), c(0, : The model should be simplified beforehand > ## leads to an error > > try(parampositions(model="whi", + param=cbind(c(1, 1, 1), c(2, 0, 0), c(1, 0, 0)))) Error in parampositions(model = "whi", param = cbind(c(1, 1, 1), c(2, : The model should be simplified beforehand > ## leads to an error > > try(parampositions(model="whi", + param=cbind(c(1, 1, 1), c(NA, 0, 0), c(1, 0, 0)))) model: whittlematern param: Error in parampositions(model = "whi", param = cbind(c(1, 1, 1), c(NA, : Model is too complex to be identified > ## leads to an error > > parampositions(model=list(list(model="exp", var=3, scale=6), "+", + list(model="whittle", var=2, scale=7, + kappa=NA))) ## again the values do not matter List of 2 $ model:List of 3 ..$ :List of 4 .. ..$ model : chr "exponential" .. ..$ var : int 1 .. ..$ kappas: NULL .. ..$ scale : int 2 ..$ op: chr "+" ..$ :List of 4 .. ..$ model : chr "whittlematern" .. ..$ var : int 3 .. ..$ kappas: int 4 .. ..$ scale : int 5 $ mean : logi NA > > > > cleanEx(); ..nameEx <- "pokeTBM" > > ### * pokeTBM > > flush(stderr()); flush(stdout()) > > ### Name: pokeTBM > ### Title: Transfer of initialisation details > ### Aliases: pokeTBM > ### Keywords: spatial > > ### ** Examples > > ###### initialisation > runif(1) [1] 0.2655087 > rs <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE) > DeleteAllRegisters() > col <- "red" > model <- "exponential" > param <- c(0, 1, 0, 10) > meth <- "TBM2" > > x2 <- seq(-50, 50, 1) > x1 <- seq(-150,150,1) > y1 <- seq(-15, 15, 1) > r1x <- range(x2) > r1y <- range(y1) > r2x <- range(x2) > r2y <- range(y1) > > ###### simulation of a random field on long thing stripe > z1 <- GaussRF(x1, y1, model=model, param=param, grid=TRUE, register=1, + method=meth) m[0]=1250, mtot=1250 using approximating circulant embedding: largest modulus of the imaginary part has been 1.813691e-14 return code 0 : Message: fine. WARNING! Positive `addradius' may have bad simulation results as consequence, difficult to detect. > get(getOption("device"))(height=1.55, width=12) > par(mar=c(2.2, 2.2, 0.1, 0.1)) > image(x1, y1, z1, col=rainbow(100)) > polygon(r1x[c(1,2,2,1)], r1y[c(1,1,2,2)], border=col) > > ###### definition of a random field on a square > InitGaussRF(x2, x2, model=model, param=param, grid=TRUE, register=2, + method=meth) m[0]=576, mtot=576 using approximating circulant embedding: largest modulus of the imaginary part has been 8.588280e-15 return code 0 : Message: fine. [1] 0 > > ###### carry over the TBM details from the stripe > if (err <- pokeTBM(1, 2)) stop(paste(err)) > > ##### simulate with the same random seed and the identical TBM details > assign(".Random.seed", rs, envir=.GlobalEnv) > z2 <- DoSimulateRF(register=2) > get(getOption("device"))(height=4.3, width=4.3) > par(mar=c(2.2, 2.2, 0.1, 0.1)) > image(x2, x2, z2, zlim=range(z1), col=rainbow(100)) > polygon(r2x[c(1,2,2,1)], r2y[c(1,1,2,2)], border=col) > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "soil" > > ### * soil > > flush(stderr()); flush(stdout()) > > ### Name: soil > ### Title: Soil data of North Bavaria, Germany > ### Aliases: soil > ### Keywords: datasets > > ### ** Examples > > > ################################################################ > ## ## > ## a geostatistical analysis that demonstrates ## > ## features of the package `RandomFields' ## > ## ## > ################################################################ > > data(soil) > names(soil) [1] "x.coord" "y.coord" "nr" "moisture" "NO3.N" "Total.N" "NH4.N" [8] "DOC" "N20N" > pts <- soil[,c(1,2)] > d <- soil$moisture > > ## define some graphical parameters first > close.screen(close.screen()) [1] FALSE > maxbin <- max(dist(pts)) / 2 > (zlim <- range(d)) [1] 4.566667 22.300000 > cn <- 100 > colour <- rainbow(cn) > par(bg="white",cex=1, cex.lab=1.4, cex.axis=1.4, mar=c(4.3,4.3,0.8,0.8)) > lu.x <- min(soil$x) > lu.y <- max(soil$y) > y <- x <- seq(min(soil$x), max(soil$x), l=121) > > ## ... and a certain appearance of the legend > my.legend <- function(lu.x, lu.y, zlim, col, cex=1) { + ## uses already the legend code of R-1.3.0 + cn <- length(col) + filler <- vector("character", length=(cn-3)/2) + legend(lu.x, lu.y, y.i=0.03, x.i=0.1, + legend=c(format(zlim[2], dig=2), filler, + format(mean(zlim), dig=2), filler, + format(zlim[1], dig=2)), + lty=1, col=rev(col),cex=cex) + } > > ## plot the data first > plot(pts, col=colour[1+(cn-1)*((d-min(d))/diff(zlim))], pch=16, + xlab="x [cm]", ylab="y [cm]", cex.axis=1.5, cex.lab=1.5) > my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1.5) > > ## empirical variogram > ev <- EmpiricalVariogram(pts, data=d, grid=FALSE, + bin=c(-1,seq(0,maxbin,l=30))) > > ## show all models, > RFparameters(Print=1, CE.force=TRUE, CE.trials=1, CE.userfft=TRUE) > by.eye <- ShowModels(x=x, y=y, emp=ev, col=colour, Zlim=zlim, + Mean=mean(d), me="ci") warning: end of action list; switching to terminal mode action=`none' > > ## fit parameters of the whittlematern model by MLE > fit <- fitvario(x=pts, data=d, model="whittle", par=rep(NA,5), + mle.m="ml", cross.m=NULL)$variogram ********+Error in chol(CC <- matrix(.C("CovarianceMatrixNatSc", distances, lc, : the leading minor of order 110 is not positive definite + > > ## plot the fitted model and the empirical variogram > plot(ev$c, ev$emp.var, ylim=c(0,11), ylab="variogram", xlab="lag") > gx <- seq(0.001, max(ev$c), l=100) > if(!is.null(by.eye)) lines(gx, Variogram(gx, model=by.eye)) > lines(gx, Variogram(gx, model=fit$ml), col=2) > lines(gx, Variogram(gx, model=fit$plain), col=3) > lines(gx, Variogram(gx, model=fit$sqrt.nr), col=4) > lines(gx, Variogram(gx, model=fit$sd.inv), col=5) > legend(120, 4, c("empirical", "by eye", "ML", "lsq", "sqrt(n) lsq", + "sd^-1 lsq"), + lty=c(-1, rep(1, 5)), pch=c(1, rep(-1, 5)), + col=c(1, 1, 2, 3, 4, 5), cex=1.4) > > ## map of expected values > k <- Kriging("O", x=x, y=y, grid=TRUE, model=fit$ml, given=pts, data=d) ............................................................................... > par(mfrow=c(1,2)) > plot(pts, col=colour[1+99*((d-min(d))/diff(zlim))], pch=16, xlab="x [cm]", ylab="y [cm]") > my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=1) > image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]") > par(bg="white") > > ## what is the probability that at no point of the > ## grid given by x and y the moisture is greater than 24 percent? > RFparameters(Print=1, CE.force=FALSE, CE.trials=3, CE.userfft=TRUE) > cs <- CondSimu("O", x=x, y=y, grid=TRUE, model=fit$ml, given=pts, + data=d, n=10) # better n=100 or n=1000 ********** > par(mfrow=c(2,3)) > image(x, y, k, col=colour, zlim=zlim, xlab="x [cm]", ylab="y [cm]") > my.legend(lu.x, lu.y, zlim=zlim, col=colour, cex=0.5) > for (i in 1:5) + image(x, y, cs[, , i], col=colour, zlim=zlim, + xlab="x [cm]", ylab="y [cm]") > mean(apply(cs<=24, 3, all)) ## about 40 percent ... [1] 0.1 > > > > graphics::par(get("par.postscript", env = .CheckExEnv)) > cleanEx(); ..nameEx <- "useraction" > > ### * useraction > > flush(stderr()); flush(stdout()) > > ### Name: useraction > ### Title: Set input behaviour > ### Aliases: useraction putactions > ### Keywords: utilities > > ### ** Examples > > sample.interaction <- function(read.r2=TRUE, type="p", pch=16) { + r1 <- Readline(prompt="first input: ") + r2 <- if (read.r2) Readline(prompt="second input: ", info="2nd input") + else " -- no input -- " + cat("Mark some points with the right mouse key, then leave with the left mouse key\n") + get(getOption("device"))() + plot(Inf, Inf, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="") + l <- Locator(100, info="locator input", type=type, pch=pch) + r3 <- Readline(prompt="third input: ", info="last input") + return(list(r1=r1, r2=r2, l=l, r3=r3)) + } > > ######################### > ## user terminal input ## > ######################### > { + useraction("start.register") + str(sample.interaction()) + str(l <- getactionlist()) + } first input: second input: ################# Mark some points with the right mouse key, then leave with the left mouse key third input: ## just replay ## List of 4 $ r1: chr "" $ r2: chr "#################" $ l : NULL $ r3: chr "## just replay ##" List of 4 $ :List of 1 ..$ : chr "" $ :List of 2 ..$ : chr "#################" ..$ info: chr "2nd input" $ :List of 1 ..$ info: chr "locator input" $ :List of 2 ..$ : chr "## just replay ##" ..$ info: chr "last input" > ################# > useraction("replay", sleep=2, wait=interactive() * 0.2 * 5, + PrintLevel=2, actionlist=l) > str(sample.interaction(type="l")) action #1 first input: action #2 (2nd input) second input: ################# Mark some points with the right mouse key, then leave with the left mouse key action #3 (locator input) action #4 (last input) third input: ## just replay ## List of 4 $ r1: chr "" $ r2: chr "#################" $ l : NULL $ r3: chr "## just replay ##" > > ############################################## > ## modify first the input list, then replay ## > ############################################## > l2 <- l[-2] > l2[[1]] <- list("some other words", info="changed input") > str(l2) List of 3 $ :List of 2 ..$ : chr "some other words" ..$ info: chr "changed input" $ :List of 1 ..$ info: chr "locator input" $ :List of 2 ..$ : chr "## just replay ##" ..$ info: chr "last input" > useraction("replay", sleep=1, wait=interactive() * 0.1, + PrintLevel=0, l2) > str(sample.interaction(read.r2=FALSE, type="o")) # input now from l2 first input: some other words Mark some points with the right mouse key, then leave with the left mouse key third input: ## just replay ## List of 4 $ r1: chr "some other words" $ r2: chr " -- no input -- " $ l : NULL $ r3: chr "## just replay ##" > Readline(prompt="new input: ", info="?!") # switch to terminal warning: end of action list; switching to terminal mode action=`none' new input: ## since end of stored list [1] "## since end of stored list" > # str(getactionlist()) # new input has not been not stored ... > > #################################################### > ## use of the two lists, l and l2, in a mixed way ## > #################################################### > useraction("replay", sleep=2, wait=interactive() * 0.2, + PrintLevel=0, actionlist=l) > Readline(prompt="first input of 1: ") first input of 1: [1] "" > dump <- getactions() > useraction("replay", sleep=0.5, wait=interactive() * 0.05, + PrintLevel=0, actionlist=l2) > Readline(prompt="first input of 2: ") first input of 2: some other words [1] "some other words" > dump2 <- putactions(dump) > Readline(prompt="second input of 1: ") second input of 1: ################# [1] "#################" > ## locator call reading from list 1: > plot(Inf, Inf, xlim=c(-1,1), ylim=c(-1,1), xlab="", ylab="") > Locator(100, info="locator input", type="p") NULL > dump <- putactions(dump2) > ##locator call reading from list 2: > Locator(100, info="locator input", type="p", pch=20, col="blue") NULL > Readline(prompt="last input of 2: ") last input of 2: ## just replay ## [1] "## just replay ##" > putactions(dump) > Readline(prompt="last input of 1: ") last input of 1: ## just replay ## [1] "## just replay ##" > > ##################################################### > ## see help("eval.parameters") for another example ## > ##################################################### > > > > > ### *