Bubble model for arbitrary areas of scales
A model that allows for arbitray areas of scale applied to an isotropic model, i.e.
C(x,y) = φ(\|x -y \| / s)
as long as s_x = s_y = s. Here, s_x is the scaling at location x,
The cross-correlations between areas of different scales are given through a modified distance d. Let z_{s} be a finite subset of R^d depending on the scale s. Let w_u be a weight for an auxiliary point u\in z_{s} with ∑_{u \in z_s} w_u = 1. Let τ_x = s_x^{-2}. Then
d^2(x, y) = \min\{τ(x), τ(y)\} \|x - y\|^2 + ∑_{ξ \in_{span(τ(x), τ(y))}} ∑_{u \in z_{ξ^{-0.5}}} w_u \|x - u\|^2 Δ ξ
Here, span(τ(x), τ(y)) is the finite set of values s^{-2} that are realized on the locations of interest and Δ ξ is the difference of two realized and ordered values of the scaling s.
RMbubble(phi, scaling, z, weight, minscale, barycentre,
var, scale, Aniso, proj)phi |
isotropic submodel |
scaling |
model that gives the non-stationary scaling s_x |
z |
matrix of the union of all z_s. The number of rows equals the dimension of the field. If not given, the locations with non-vanishing gradient are taken. |
weight |
vector of weights w
whose length equals the number of columns of |
minscale |
vector for partioning z into classes z_s.
Its length equals the number of columns of |
barycentre |
logical. If The argument has no effect when Default: |
var,scale,Aniso,proj |
optional arguments; same meaning for any
|
minscale gives the minimal scale s value above which
the corresponding points z define the set z_s.
The validity of the set z_s ends with the next lower value
given.
Let minscale = (10, 10, 10, 7, 7, 7, 0.5). Then for some
d-dimensional vectors z_1,…, z_7 we have
z_s = \{ z_1, z_2, z_3 \}, s ≥ 10
z_s = \{ z_4, z_5, z_5 \}, 7 ≥ s < 10
z_s = \{ z_7 \}, s ≥ 0.5
Note that, in this case, all realized scaling values must be ≥ 0.5. Note further, that the weights for the subset must sum up to one, i.e.
w_1+w_2 +w_3=w_4 + w_5 + w_6 = w_7 = 1.
This model is defined only for grids.
Martin Schlather, schlather@math.uni-mannheim.de, https://www.wim.uni-mannheim.de/schlather/
Bonat, W.H. , Ribeiro, P. Jr. and Schlather, M. (2019) Modelling non-stationarity in scale. In preparation.
RFoptions(seed=0) ## *ANY* simulation will have the random seed 0; set
## RFoptions(seed=NA) to make them all random again
x <- seq(0,1, if (interactive()) 0.02 else 0.5)
d <- sqrt(rowSums(as.matrix(expand.grid(x-0.5, x-0.5))^2))
d <- matrix(d < 0.25, nc=length(x))
image(d)
scale <- RMcovariate(data=as.double(d) * 2 + 0.5, raw=TRUE)
## two models:
## the frist uses the standard approach for determining the
## reference point z, which is based on gradients
## the second takes the centre of the ball
model1 <- RMbubble(RMexp(), scaling=scale)
model2 <- RMbubble(RMexp(), scaling=scale, z=c(0.5, 0.5))
model3 <- RMbubble(RMexp(), scaling=scale, barycentre=TRUE) # approx. of model2
## model2 has slightly higher correlations than model1:
C1 <- RFcovmatrix(model1, x, x)
C2 <- RFcovmatrix(model2, x, x)
C3 <- RFcovmatrix(model3, x, x)
print(range(C2 - C1))
dev.new(); hist(C2 - C1)
print(range(C3 - C2)) # only small differences to C2
print(mean(C3 - C2))
dev.new(); hist(C3 - C2)
plot(z1 <- RFsimulate(model1, x, x))
plot(z2 <- RFsimulate(model2, x, x))
plot(z3 <- RFsimulate(model3, x, x)) # only tiny differences to z2
## in the following we compare the standard bubble model with
## the models RMblend, RMscale and RMS (so, model2 above
## performs even better)
biwm <- RMbiwm(nudiag=c(0.5, 0.5), nured=1, rhored=1, cdiag=c(1, 1),
s=c(0.5, 2.5, 0.5))
blend <- RMblend(multi=biwm, blend=RMcovariate(data = as.double(d), raw=TRUE))
plot(zblend <- RFsimulate(blend, x, x)) ## takes a while ...
Cblend <- RFcovmatrix(blend, x, x)
Mscale <- RMscale(RMexp(), scaling = scale, penalty=RMid() / 2)
plot(zscale <- RFsimulate(Mscale, x, x))
Cscale <- RFcovmatrix(Mscale, x, x)
Mscale2 <- RMscale(RMexp(), scaling = scale, penalty=RMid() / 20000)
plot(zscale2 <- RFsimulate(Mscale2, x, x))
Cscale2 <- RFcovmatrix(Mscale2, x, x)
S <- RMexp(scale = scale)
plot(zS <- RFsimulate(S, x, x))
CS <- RFcovmatrix(S, x, x)
print(range(C1 - CS))
print(range(C1 - Cscale))
print(range(C1 - Cscale2))
print(range(C1 - Cblend))
dev.new(); hist(C1-CS) ## C1 is better
dev.new(); hist(C1-Cscale) ## C1 is better
dev.new(); hist(C1-Cscale2) ## both are equally good. Maybe C1 slightly better
dev.new(); hist(C1-Cblend) ## C1 is betterPlease choose more modern alternatives, such as Google Chrome or Mozilla Firefox.