Estimating Fractal Dimensions

Description

Implements various methods for estimating fractal dimension of time series and 2-dimensional data.

Details

The package provides tools for estimating fractal dimension of one- or two-dimensional data, using methods described in Gneiting et al. (2010). The user can take an advantage of the available sliding window technique in which a window of a given size is slided along the data and an estimate is obtained for each position.

The main function is fd.estimate which can be used for one dimensional time series, as well as for two dimensional data. It computes one estimate for each method and each sliding window. It is a wrapper for lower level functions for computing just one estimate on the given data, see fd.estim.method for details.

Author(s)

Hana Sevcikova <hanas@uw.edu>, Tilmann Gneiting <tilmann@stat.washington.edu>, Don Percival <dbp@apl.washington.edu>

Maintainer: Hana Sevcikova <hanas@uw.edu>

References

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. <doi:10.1214/11-STS370> (Version as technical report available at https://stat.uw.edu/sites/default/files/files/reports/2010/tr577.pdf)

See Also

fd.estimate, fd.estim.method

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

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:

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

Estimating Fractal Dimensions of Time Series and Two-dimensional Data

Description

The functions compute a set of fractal dimensions D for time series and two-dimensional data via various methods using a sliding window technique. There is one D computed for each method and for each sliding window of a given size that is moved along the data.

Usage

## S3 method for class 'numeric'
 fd.estimate(data, methods = "madogram", window.size = length(data), 
    step.size = window.size, trim = TRUE, keep.data = FALSE, 
    keep.loglog = FALSE, parallel = FALSE, nr.nodes = NULL, 
    debuglevel = 0, ...)
	
## S3 method for class 'matrix'
 fd.estimate(data, methods = "transect.var", window.size = ncol(data),
    step.size = window.size, trim = TRUE, keep.data = FALSE, 
    keep.loglog = FALSE, parallel = FALSE, nr.nodes = NULL, 
    debuglevel = 0, ...)
          

Arguments

Details

In case of one-dimensional time series, the function initiates a sliding window of the given size at the beginning of the time series. The window is moved along the data by the given step size. If parallel is TRUE, computation on each window happens in parallel. In the two-dimensional case, the window is initiated in the top left corner of the data matrix and moved horizontally by the given step size, as well as vertically by the same step size. If the process is running in parallel, processing each row is done in parallel. In both cases, in each iteration estimates of fractal dimension for data within the sliding window are computed using the given estimation methods.

Note that the estimation results are NA for any sliding window that contains NA values.

Arguments that are to be passed to specific methods can be given either directly, if they applies to all given methods. Or, they can be given as a list via the methods argument: There is one list per method that must contain the entry “name” being the method name. Remaining entries in the list corespond to one argument each (see Example below).

Value

An object of class FractalDim which consists of the following components:

See Also

fd.estim.method, fd.get.available.methods, FDloglog, fd.get

Examples

## Not run: 
library(RandomFields)
n <- 10000
# generate a time series
rf <- GaussRF(x = c(0, 1, 1/n), model = "stable", 
       grid = TRUE, gridtriple = TRUE,
       param = c(mean=0, variance=1, nugget=0, scale=100, kappa=1))

# Plots for two sliding windows of each of the four methods below.
# Argument nlags is common to all methods;
# the 'variation' method has in addition argument p.index
par(mfrow=c(2,4)) # one row per window
fd <- fd.estimate(rf, 
       methods = list(list(name="variation", p.index=0.5), 
       					"variogram", "hallwood", "boxcount"),
       window.size = 5000, step.size = 5000, plot.loglog = TRUE, nlags = 10)
    
# 2d random fields
n <- 200
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(2,2))
# plots for 4 sliding windows (2 horizontal, 2 vertical)
fd2d <- fd.estimate(rf2d, methods="filter1",
         window.size = 100, step.size=100, plot.loglog = TRUE)

## End(Not run)

Access Method for Objects of Class FractalDim

Description

For given method it returns the corresponding estimates.

Usage

fd.get(fractaldim, method)

Arguments

Value

A FractalDim object. The original fractaldim arrays fd and scale are reduced in the last dimension into only one method, namely the given method.

See Also

fd.estimate, fd.get.available.methods

Examples

## Not run: 
library(RandomFields)
x <- seq(0, 10000)
# generate a random field
truealpha <- 1.5
rf <- GaussRF(x = x, model = "stable", grid = TRUE,
    param = c(mean=0, variance=1, nugget=0, scale=100, 
        alpha=truealpha))

#compute fractal dimension using various methods
methods <- c("madogram", "variogram", "hallwood", "boxcount",
                       "periodogram","dctII", "wavelet")
fdts <- fd.estimate (rf, methods = methods, window.size = 500, 
    step.size = 100, nlags = 10, trim = FALSE, debuglevel = 3)

# plot the variation
cols <- rainbow(length(methods))
plot(ts(fd.get (fdts, methods[1])$fd),ylim=c(min(fdts$fd), max(fdts$fd)),
             ylab="fd", col=cols[1])
for (imeth in 2:length(methods)) 
    lines(ts(fd.get (fdts, methods[imeth])$fd), col=cols[imeth])
legend('topleft', legend=methods, col=cols, lwd=1)
abline(h=2-truealpha/2)

## End(Not run)

Available Estimation Methods

Description

The function returns a list of estimation methods that can be used in fd.estimate and other functions of the package.

Usage

fd.get.available.methods(dim = 1)

Arguments

Value

A list of the available methods. Their order number in the list corresponds to the internal codes of the methods.

See Also

fd.estimate, fd.get

Internal Functions

Description

Internal functions.

Details

Functions for internal use only.

Obtaining Regression Object

Description

Obtaining and summarizing result of a linear regression of the log-log plot, object of class FDloglog.

Usage

get.rawFD.from.regression(x, y, leaveout = 0)

## S3 method for class 'FDloglog'
summary(object, ...)

Arguments

Value

Function get.rawFD.from.regression returns an object of class FDloglog with the following components:

Author(s)

Hana Sevcikova

See Also

fd.estimate, fd.estim.method

Summary for an Object of Class FractalDim

Description

The function prints summary of estimates in a FractalDim object.

Usage

## S3 method for class 'FractalDim'
 summary(object, ...)

Arguments

Details

The function prints information about the grid on which the estimates were obtained. For each method it shows the mean and standard deviation of data in each of the two components (fd and scale).

See Also

fd.estimate