The functions estimate a fractal dimension of the given data. Each function uses a different method. Functions for boxcount, hallwood, variogram, madogram, rodogram, variation, incr1, genton, periodogram, wavelet and dctII methods are to be used on one-dimensional time series. The remaining functions (transect, isotropic, squareincr, and filter1) are to be used on two-dimensional data.
fd.estim.boxcount (data, plot.loglog = FALSE, nlags = "auto",
shift.up=TRUE, plot.allpoints = FALSE, legend.type = 's',
…, debuglevel = 0)
fd.estim.hallwood (data, plot.loglog = FALSE, nlags = "auto",
plot.allpoints = FALSE, legend.type = 's', …, debuglevel = 0)
fd.estim.variogram (data, …)
fd.estim.madogram (data, …)
fd.estim.rodogram (data, …)
fd.estim.variation (data, p.index = 1, …)
fd.estim.incr1(data, p.index=2, …)
fd.estim.genton (data, …)
fd.estim.periodogram (data, plot.loglog = FALSE, nlags = "auto", …)
fd.estim.wavelet (data, plot.loglog=FALSE, plot.allpoints = FALSE,
filter = "haar", J1 = max(1,floor(log2(length(data))/3-1)),
J0 = floor(log2(length(data))), legend.type = 's',
…, debuglevel = 0)
fd.estim.dctII (data, plot.loglog = FALSE, nlags = "auto", …)
fd.estim.transect.var (data, p.index = 2, …)
fd.estim.transect.incr1 (data, p.index = 2, …)
fd.estim.isotropic (data, p.index = 2, direction = 'hvd+d-',
plot.loglog = FALSE, nlags = "auto", plot.allpoints = FALSE,
legend.type = 's', …, debuglevel=0)
fd.estim.squareincr (data, p.index = 2,
plot.loglog = FALSE, nlags = "auto", plot.allpoints = FALSE,
legend.type = 's', …, debuglevel=0)
fd.estim.filter1 (data, p.index = 2, direction = 'hvd+d-',
plot.loglog = FALSE, nlags = "auto", plot.allpoints = FALSE,
legend.type = 's', …, debuglevel=0)For the first eleven functions data is a one-dimensional vector. For the last five functions data is a matrix.
Parameter \(p\) of the variation method (see below).
For the 2d estimators, this argument specifies the direction of the estimation (see details below). It can be any combination of the characters ‘h’ (horizontal), ‘v’ (vertical), ‘d+’ (diagonal with positive gradient), and ‘d-’ (diagonal with negative gradient). These characters should be combined into one string.
Logical value determining if the underlying log-log plots should be plotted.
Number of lags to be used in the
estimation. Possible values are "auto", "all" or a
single number. If nlags = "auto", each method sets the number of lags
to the theoretically "best" value for that method. "all" means
that all lags are included in the estimation.
For each interval on the horizontal axis, it moves the boxes vertically up to the smallest data point of that interval. If it is FALSE, all boxes are on a regular grid.
Logical. If FALSE, only points that were considered in
the regression are shown. Otherwise, all points of the log-log
plot are shown in the graph and those considered in the regression are marked by filled circles.
This argument is only used if
plot.loglog = TRUE. Note that setting this argument to TRUE might (depending on the method)
considerably increase the computation run-time.
Argument passed to the modwt function of the wavelets package.
Parameters of the wavelet method controlling the number of frequencies used in the estimation.
One of the characters 'f', 's', or 'n'. It controlls the amount of information in the legend of the log-log plot. If it is 'f' (full), values of fd and scale, including the raw values of the corresponding slope and intercept are shown. If it is 's' (short), only fd is shown. Value of 'n' (None) causes no legend being plotted. The argument is only used if plot.loglog = TRUE.
Arguments passed to the plotting function if plot.loglog = TRUE. For some functions, … contain additional arguments, see Details.
Controls the amount of debugging messages. The functions produce messages on level 5.
Each function returns an object of class FractalDim with elements:
Here it is always 1.
Single value, namely the estimated fractal dimension and scale, respectively.
Method name and code used for the estimation.
Size of the data.
Dimension of the data used for the estimation. It is either one or two.
Object of class FDloglog used for the estimation.
The methodology of these functions is based on the theory described in Gneiting et al (2010). Please refer to this paper for notation. Here we give only a few comments about the implementation.
The function fd.estim.boxcount determines the smallest possible value of \(m\)
for which \(n\leq 2^m\) is a power of 2. Only data points \(x_1,\dots,x_{n_{eff}}\) are considered for the estimation,
where \(n_{eff} = 2^{m-1}+1\). The value of \(K\) can be given by the user through the argument nlags. If nlags = "auto", box sizes \(\epsilon_k\) for \(k=j,j+1,\dots,m-2\) are considered,
where for all \(i<j\) is \(N(\epsilon_i) > \frac{n_{eff}}{5}\), i.e.
the two largest box sizes and very small boxes are eliminated (corresponds to the
Liebovitch and Toth modification).
If shift.up=TRUE, the algorithm shifts each vertical column of boxes up to the smallest data value in that column.
\(N(\epsilon_k)\) for a particular box is increased if either a data point is contained in the box, or if a line connecting two neighboring data points crosses the box.
This estimator is a version of box-count that
instead of number of boxes considers the area of boxes that cover
the underlying curve.
Hall and Wood (1993) recommend
the use of \(L = 2\) which the function fd.estim.hallwood uses
if the arguments nlags = "auto".
The p.index argument of fd.estim.variation and fd.estim.incr1 is the power index \(p\). The madogram, variogram, and rodogram, respectively, correspond to the Variation estimator with \(p\) equals 1, 2, and 1/2, respectively. The Incr1 estimator is like Variation but based on second order differences.
Any argument that can be passed to fd.estim.hallwood can be passed here as well. In addition,
as in the Hall-Wood case, \(L\) is set to 2 for these estimators, if nlags = "auto".
This is a highly robust variogram estimator as proposed by Genton (1998). Given \(U_i(d) = X_{i/n} - X_{(i-d)/n}\), define $$ \hat{V}(d) = \left[2.2191\{|U_i(d) - U_j(d)|; i<j \}_{(k)}\right]^2 , \quad \mbox{where} \;\; k={\lfloor (n-d)/2\rfloor +1 \choose 2}\,.$$
Thus, the estimator is derived from the \(k\)-th quantile of the \(U_i(d)\) values. The \(\hat{D}_k\) estimator is derived from the log-log plot of \(\log(d)\) against \(\log(\hat{V}(d))\). The implementation uses the qn function of the pcaPP package to compute \(\hat{V}(d)\).
Here again, the number of lags is set to 2 if nlags = "auto" and any arguments of the fd.estim.hallwood are accepted here as well.
The method is implemented as proposed by Chan et al. (1995) with notation from Gneiting et al (2010).
As Chan et al. (1995) recommend, we use \(L=\lfloor \min(m/2, n^{2/3})\rfloor\) if nlags = "auto". Any arguments of the fd.estim.hallwood are also accepted here.
This method uses \(J_0\) vectors of wavelet coefficients which are obtained using the function modwt of the wavelets package. The choice of J0 and J1 determine the number of frequencies used in the estimation.
If nlags = "auto", we use \(L=\lfloor \min(2m, 4n^{2/3})\rfloor\). Any arguments of the fd.estim.hallwood are also accepted here.
The two-dimensional estimators are all based on the Variation method with the power index \(p\) (argument p.index) with the following alternatives:
For every given direction, a variation estimate (or a variant that uses second differences) is found in each row (for horizontal direction) and/or column (for vertical direction). The resulting estimate is the median over the set of estimates. In the function fd.estim.transect.var the line transect estimates are based on first differences; In the function fd.estim.transect.incr1 they are based on second differences.
This method does not support the feature of creating a log-log plot, since there are many log-log regressions from which the results are derived. The methods also accept arguments direction, nlags and debuglevel.
Davies and Hall (1999) on page 12 define the isotropic empirical variogram. This is here implemented more generally using the variation estimator. If nlags = "auto", the number of lags is set to either 3 if diagonal direction is used together with either horizontal or vertical direction or both. If only horizontal or/and vertical direction is used, the number of lags is set to 2.
We use the square-increment estimator proposed in eqs. (4.2) through (4.7) of Chan and Wood (2000). Note that this method is equivalent to the Filter 3 approach of Zhu and Stein (2002) which is the way it is implemented in the package. The automatic setting of number of lags is done as for the Isotropic method.
Here, the Filter 1 approach of Zhu and Stein (2002) is implemented. Again, the automatic setting of number of lags is done as for the Isotropic method.
For all methods (but Transect), if the argument plot.loglog is TRUE, a graph with the log-log plot is shown,
including the fitted regression line. Only points included in the regression are plotted, unless the argument plot.allpoints is set to TRUE. In such a case, points used for fitting the regression line are marked by filled circles.
For using multiple estimation methods via one function see fd.estimate.
Chan, G., Hall, P., Poskitt, D. (1995) Periodogram-Based Estimators of Fractal Properties. Annals of Statistics 23 (5), 1684--1711.
Chan, G., Wood, A. (2000) Increment-based estimators of fractal dimension for two-dimensional surface data. Statistica Sinica 10, 343--376.
Davies, S., Hall, P. (1999) Fractal analysis of surface roughness by using spatial data. Journal of the Royal Statistical Society Series B 61, 3--37.
Genton, M. G. (1998) Highly robust variogram estimation. Mathematical Geology 30, 213--221.
Gneiting, T., Sevcikova, H. and Percival, D. B. (2012). Estimators of fractal dimension: Assessing the smoothness of time series and spatial data. Statistical Science, 27(2), 247-277. (Version as technical report available at https://stat.uw.edu/sites/default/files/files/reports/2010/tr577.pdf)
Hall, P., Wood, A. (1993) On the Performance of Box-Counting Estimators of Fractal Dimension. Biometrika 80 (1), 246--252.
Zhu, Z., Stein, M. (2002) Parameter estimation for fractional Brownian surfaces. Statistica Sinica 12, 863--883.
# NOT RUN {
if (requireNamespace("RandomFields", quietly = TRUE)) withAutoprint({
library(RandomFields)
# 1d time series
n <- 256
rf <- GaussRF(x = c(0,1, 1/n), model = "stable",
grid = TRUE, gridtriple = TRUE,
param = c(mean=0, variance=1, nugget=0, scale=1, kappa=1))
par(mfrow=c(4,2))
fd.estim.variogram (rf, nlags = 20, plot.loglog = TRUE)
fd.estim.variation (rf, nlags = 20, plot.loglog = TRUE)
fd.estim.variogram (rf, nlags = 3, plot.loglog = TRUE,
plot.allpoints = TRUE)
fd.estim.variation (rf, plot.loglog = TRUE, plot.allpoints = TRUE)
fd.estim.hallwood (rf, nlags = 10, plot.loglog = TRUE)
fd.estim.boxcount (rf, nlags = "all", plot.loglog = TRUE,
plot.allpoints = TRUE)
fd.estim.periodogram (rf, plot.loglog = TRUE)
fd.estim.dctII (rf, plot.loglog = TRUE)
# 2d random fields
n <- 128
rf2d <- GaussRF(x = c(0,1, 1/n), y = c(0,1, 1/n), model = "stable",
grid = TRUE, gridtriple = TRUE,
param = c(mean=0, variance=1, nugget=0, scale=1, kappa=1))
par(mfrow=c(1,3))
fd.estim.isotropic (rf2d, p.index = 1, direction='hv',
plot.loglog = TRUE, plot.allpoints = TRUE)
fd.estim.squareincr (rf2d, p.index = 1, plot.loglog = TRUE, plot.allpoints = TRUE)
fd.estim.filter1 (rf2d, p.index = 1, plot.loglog = TRUE, plot.allpoints = TRUE)
})
# }
Run the code above in your browser using DataLab