fd.estim.method {fractaldim}

Estimation of Fractal Dimension via Specific Methods

Description

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.

Usage

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)

Arguments

data

For the first eleven functions data is a one-dimensional vector. For the last five functions data is a matrix.

p.index

Parameter p of the variation method (see below).

direction

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.

plot.loglog

Logical value determining if the underlying log-log plots should be plotted.

nlags

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.

shift.up

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.

plot.allpoints

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.

filter

Argument passed to the modwt function of the wavelets package.

J0, J1

Parameters of the wavelet method controlling the number of frequencies used in the estimation.

legend.type

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.

debuglevel

Controls the amount of debugging messages. The functions produce messages on level 5.

Details

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.

Box-count estimator:

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.

Hall-Wood estimator

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".

Variation, Variogram, Madogram, Rodogram, and Incr1 estimators:

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".

Genton robust estimator:

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.

Periodogram estimator:

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.

Wavelet estimator:

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.

DCT-II estimator:

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:

Transect

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.

Isotropic

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.

Square-increment

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.

Filter 1

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.

Value

Each function returns an object of class FractalDim with elements:

dim

Here it is always 1.

fd, scale

Single value, namely the estimated fractal dimension and scale, respectively.

methods, methods.coding

Method name and code used for the estimation.

window.size, step.size

Size of the data.

data.dim

Dimension of the data used for the estimation. It is either one or two.

loglog

Object of class FDloglog used for the estimation.

Note

Function fd.estimate can be used as a wrapper for these functions.

Author(s)

Hana Sevcikova, Don Percival, Tilmann Gneiting

References

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.

See Also

fd.estimate

Examples

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)
})