| Title: | Continuous-Lag Spatial Markov Chains |
|---|---|
| Description: | A set of functions is provided for 1) the stratum lengths analysis along a chosen direction, 2) fast estimation of continuous lag spatial Markov chains model parameters and probability computing (also for large data sets), 3) transition probability maps and transiograms drawing, 4) simulation methods for categorical random fields. More details on the methodology are discussed in Sartore (2013) <doi:10.32614/RJ-2013-022> and Sartore et al. (2016) <doi:10.1016/j.cageo.2016.06.001>. |
| Authors: | Luca Sartore [aut, cre]
|
| Maintainer: | Luca Sartore <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.3.15 |
| Built: | 2024-09-01 04:01:16 UTC |
| Source: | https://github.com/cran/spMC |
The main goal of this package is to provide a set of functions for
the stratum lengths analysis along a chosen direction,
fast estimation of continuous lag spatial Markov chains model parameters and probability computing (also for large data sets),
transition probability maps and transiograms drawing,
simulation methods for categorical random fields.
| Package: | spMC |
| Type: | Package |
| Version: | 0.3.15 |
| Date: | 2023-04-30 |
| License: | GPL (>= 2) |
| LazyLoad: | yes |
Several functions are available for the stratum lengths analysis, in particular they compute the stratum lengths for each stratum category, they compute the empirical distributions and many other tools for a graphical analysis.
Usually, the basic inputs for the most of the functions are a vector of categorical data and their location coordinates. They are used to estimate empirical transition probabilities (transiogram), to estimate model parameters (tpfit for one-dimensional Markov chains or multi_tpfit for multidimensional Markov chains). Once parameters are estimated, it's possible to compute theoretical transition probabilities by the use of the function predict.tpfit for one-dimensional Markov chains and predict.multi_tpfit for multidimensional ones.
The function plot.transiogram allows to plot one-dimensional transiograms, while image.multi_tpfit permit to draw transition probability maps. A powerful tool to explore graphically the anisotropy of such process is given by the functions pemt and image.pemt, which let the user to draw "quasi-empirical" transition probability maps.
Simulation methods are based on Indicator Kriging (sim_ik), Indicator Cokriging (sim_ck), Fixed or Random Path algorithms (sim_path) and Multinomial Categorical Simulation technique (sim_mcs).
Luca Sartore
Maintainer: Luca Sartore [email protected]
Allard, D., D'Or, D., Froidevaux, R. (2011) An efficient maximum entropy approach for categorical variable prediction. European Journal of Soil Science, 62(3), 381-393.
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Dynkin, E. B. (1961) Theory of Markov Processes. Englewood Cliffs, N.J.: Prentice-Hall, Inc.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Li, W. (2007) A Fixed-Path Markov Chain Algorithm for Conditional Simulation of Discrete Spatial Variables. Mathematical Geology, 39(2), 159-176.
Li, W. (2007) Markov Chain Random Fields for Estimation of Categorical Variables. Mathematical Geology, 39(June), 321-335.
Li, W. (2007) Transiograms for Characterizing Spatial Variability of Soil Classes. Soil Science Society of America Journal, 71(3), 881-893.
Pickard, D. K. (1980) Unilateral Markov Fields. Advances in Applied Probability, 12(3), 655-671.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
Sartore, L. (2013). spMC: Modelling Spatial Random Fields with Continuous Lag Markov Chains. The R Journal, 5(2), 16-28.
Sartore, L., Fabbri, P. and Gaetan, C. (2016). spMC: an R-package for 3D lithological reconstructions based on spatial Markov chains. Computers & Geosciences, 94(September), 40-47.
Weise, T. (2009) Global Optimization Algorithms - Theory and Application. https://archive.org/details/Thomas_Weise__Global_Optimization_Algorithms_Theory_and_Application.
The data set refers to a sampled area which is located in the province of Venice. Its sample units report the geographical position of the perforation, the depth, the ground permeability and other two categorical variables which denote the soil composition.
data(ACM)data(ACM)
A data frame with 2321 observations on the following 6 variables.
Xa numeric vector (longitude)
Ya numeric vector (latitude)
Za numeric vector (depth)
MAT5a factor with levels Clay, Gravel, Mix of Sand and Clay, Mix of Sand and Gravel and Sand
MAT3a factor with levels Clay, Gravel and Sand
PERMa logical vector (symmetric dichotomous variable)
Fabbri, P. (2010) Professor at the Geosciences Department of the University of Padua.
[email protected]
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) str(ACM) summary(ACM)data(ACM) str(ACM) summary(ACM)
Produce box-and-whisker plots of the stratum lengths.
## S3 method for class 'lengths' boxplot(x, ..., log = FALSE, zeros.rm = TRUE)## S3 method for class 'lengths' boxplot(x, ..., log = FALSE, zeros.rm = TRUE)
x |
an object of the class |
... |
other arguments to pass to the function |
log |
a logical value. If |
zeros.rm |
a logical value. If |
The box-and-whisker plots give some information about the distribution of the stratum lengths for the observed categories along a given direction.
An image is produced on the current graphics device. The function returns a list with the following components:
stats |
a matrix containing the values used to plot the box-and-whisker plots. |
n |
a vector with the number of observations for each category. |
conf |
a matrix containing further values to draw the lower and upper extremes of the notch. |
out |
a vectors with the values of the outlier points. |
group |
a vector whose elements indicate to which category the outlier belongs. |
names |
a character vector with the names of each category. |
Luca Sartore [email protected]
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Make the boxplot of the object gl boxplot(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Make the boxplot of the object gl boxplot(gl)
The function draws the -D sections contour plots of a multi-directional transiogram computed without any ellipsoidal interpolation and superpose the contour lines of the theoretical transition probabilities.
## S3 method for class 'pemt' contour(x, nlevels = 10, col = c("black", "blue"), main, mar, ask = TRUE, ...)## S3 method for class 'pemt' contour(x, nlevels = 10, col = c("black", "blue"), main, mar, ask = TRUE, ...)
x |
an object of class |
nlevels |
the number of levels to pass to the function |
col |
a vector of two colors to pass to the function |
main |
the main title (on top) whose font and size are fixed. |
mar |
a scalar or a numerical vector of the form |
ask |
a logical value; if |
... |
other arguments to pass to the function |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. The probability is computed for any lag vector through
where entries of are not ellipsoidally interpolated, but they are estimated for the direction specified by the vector .
The exponential matrix is evaluated by the scaling and squaring algorithm.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
image.pemt, contour, plot.transiogram
data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "avg") # Contour plots of 2-D sections of # multi-directional transiograms contour(psEmpTr, mar = .7)data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "avg") # Contour plots of 2-D sections of # multi-directional transiograms contour(psEmpTr, mar = .7)
The function estimates the empirical conditional density of the stratum lengths given the category.
## S3 method for class 'lengths' density(x, ..., log = FALSE, zeros.rm = TRUE)## S3 method for class 'lengths' density(x, ..., log = FALSE, zeros.rm = TRUE)
x |
an object of the class |
... |
other arguments to pass to the function |
log |
a logical value. If |
zeros.rm |
a logical value. If |
The function estimates the empirical density of the stratum lengths for each category by the use of the kernel methodology.
An object of class density.lengths is returned. It contains objects of class density, the given direction of the stratum lengths and a logical value which points out if the density is computed for the logarithm of stratum lengths.
Luca Sartore [email protected]
Simonoff, J. S. (1996) Smoothing Methods in Statistics. Springer-Verlag.
getlen, density.default, plot.density.lengths, print.density.lengths
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum lengths dgl <- density(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum lengths dgl <- density(gl)
The function estimates the embedded transition probabilities matrix for a -D spatial embedded Markov chain.
embed_MC(data, coords, loc.id, direction)embed_MC(data, coords, loc.id, direction)
data |
a categorical data vector of length |
coords |
an |
loc.id |
a vector of |
direction |
a |
An embedded Markov chain is probabilistic model which defines the transition probabilities between embedded occurrences.
The resulting matrix is given by normalizing a transition count matrix, which doesn't depend on the length of embedded occurrences. Self-transitions of embedded occurrences are not observable, so diagonal entries are set to be NA.
It's also possible to calculate the transition probabilities matrix for several directions in a -D space through arguments direction and loc.id. If the user has no previous knowledge about loc.id, the function which_lines provides a method to compute the right values.
A transition probability matrix, where denotes the number of observed categories. Another matrix with the counts of transitions is attached as an attribute.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Dynkin, E. B. (1961) Theory of Markov Processes. Englewood Cliffs, N.J.: Prentice-Hall, Inc.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
which_lines, predict.tpfit, predict.multi_tpfit
data(ACM) direction <- c(0, 0, 1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction, pi/8) # Estimate the embedded transition probabilities # matrix for the categorical variable MAT5 embed_MC(ACM$MAT5, ACM[, 1:3], loc.id, direction) # Estimate the embedded transition probabilities # matrix for the categorical variable MAT3 embed_MC(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Estimate the embedded transition probabilities # matrix for the categorical variable PERM embed_MC(ACM$PERM, ACM[, 1:3], loc.id, direction)data(ACM) direction <- c(0, 0, 1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction, pi/8) # Estimate the embedded transition probabilities # matrix for the categorical variable MAT5 embed_MC(ACM$MAT5, ACM[, 1:3], loc.id, direction) # Estimate the embedded transition probabilities # matrix for the categorical variable MAT3 embed_MC(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Estimate the embedded transition probabilities # matrix for the categorical variable PERM embed_MC(ACM$PERM, ACM[, 1:3], loc.id, direction)
The function estimates the stratum lengths for a -D spatial embedded Markov chain for a specified direction .
getlen(data, coords, loc.id, direction, zero.allowed = FALSE)getlen(data, coords, loc.id, direction, zero.allowed = FALSE)
data |
a categorical data vector of length |
coords |
an |
loc.id |
a vector of |
direction |
a |
zero.allowed |
a logical value which allows to return zero stratum lengths. It is |
Stratum lengths are the lengths occupied by the same -th category along lines in the direction .
A list containing the following components:
length |
a numerical vector with the stratum lengths along the given direction. |
categories |
a vector with the stratum categories. |
maxcens |
a vector with the maxima estimated censored lengths for each stratum. |
directions |
a |
zeros |
a logical values which denotes the possible presence of zero lengths. |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT5, ACM[, 1:3], loc.id, direction)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT5, ACM[, 1:3], loc.id, direction)
The function compute the histograms of the stratum lengths for each category. If plot = TRUE, the resulting object of class hist.lengths is plotted before it is returned.
## S3 method for class 'lengths' hist(x, ..., log = FALSE, zeros.rm = TRUE)## S3 method for class 'lengths' hist(x, ..., log = FALSE, zeros.rm = TRUE)
x |
an object of the class |
... |
further arguments to pass to the function |
log |
a logical value. If |
zeros.rm |
a logical value. If |
If plot = TRUE, an image is produced on the current graphics device. The function returns an object of class hist.lengths. It contains class histogram objects, the given direction of the stratum lengths and a logical value which points out if histograms are computed for the logarithm of stratum lengths.
Luca Sartore [email protected]
getlen, hist, density.lengths, plot.density.lengths
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Plot the histograms hist(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Plot the histograms hist(gl)
The function plots -D sections of a predicted multidimensional transiograms computed through ellipsoidal interpolation.
## S3 method for class 'multi_tpfit' image(x, mpoints, which.dire, max.dist, main, mar, ask = TRUE, ..., nlevels = 10, contour = TRUE)## S3 method for class 'multi_tpfit' image(x, mpoints, which.dire, max.dist, main, mar, ask = TRUE, ..., nlevels = 10, contour = TRUE)
x |
an object of the class |
mpoints |
the number of points per axes. It controls the accuracy of images to plot. |
which.dire |
a vector with two chosen axial directions. If omitted, all |
max.dist |
a scalar or a vector of maximum length for the chosen axial directions. |
main |
the main title (on top) whose font and size are fixed. |
mar |
a scalar or a numerical vector of the form |
ask |
a logical value; if |
... |
other arguments to pass to the function |
nlevels |
the number of levels to pass to the function |
contour |
logical. If |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. It is computed for any lag vector through
where entries of are ellipsoidally interpolated (see multi_tpfit).
The exponential matrix is evaluated by the scaling and squaring algorithm.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
multi_tpfit, pemt, image.pemt, image, plot.transiogram
data(ACM) # Estimate model parameter x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Set short names for categories 3 and 4 names(x$prop)[3:4] <- c("Clay and Sand", "Gravel and Sand") # Plot 2-D theoretical sections of # a multidimensional transiogram image(x, 40, max.dist=c(200,200,20), which.dire=2:3, mar = .7, col=rev(heat.colors(500)), breaks=0:500/500, nlevels = 5)data(ACM) # Estimate model parameter x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Set short names for categories 3 and 4 names(x$prop)[3:4] <- c("Clay and Sand", "Gravel and Sand") # Plot 2-D theoretical sections of # a multidimensional transiogram image(x, 40, max.dist=c(200,200,20), which.dire=2:3, mar = .7, col=rev(heat.colors(500)), breaks=0:500/500, nlevels = 5)
The function plots -D sections of a multidirectional transiogram computed without any ellipsoidal interpolation.
## S3 method for class 'pemt' image(x, main, mar, ask = TRUE, ..., nlevels = 10, contour = TRUE)## S3 method for class 'pemt' image(x, main, mar, ask = TRUE, ..., nlevels = 10, contour = TRUE)
x |
an object of class |
main |
the main title (on top) whose font and size are fixed. |
mar |
a scalar or a numerical vector of the form |
ask |
a logical value; if |
... |
other arguments to pass to the function |
nlevels |
the number of levels to pass to the function |
contour |
logical. If |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. The probability is computed for any lag vector through
where entries of are not ellipsoidally interpolated, but they are estimated for the direction specified by the vector .
The exponential matrix is evaluated by the scaling and squaring algorithm.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
image.multi_tpfit, image, plot.transiogram
data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "mlk") # Plot 2-D sections of # a multi-directional transiogram image(psEmpTr, col = rev(heat.colors(500)), breaks = 0:500 / 500, mar = .7, contour = FALSE)data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "mlk") # Plot 2-D sections of # a multi-directional transiogram image(psEmpTr, col = rev(heat.colors(500)), breaks = 0:500 / 500, mar = .7, contour = FALSE)
Function to test if an object is of the class lengths.
is.lengths(object)is.lengths(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a lengths object.
A logical value.
Luca Sartore [email protected]
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Test the object gl is.lengths(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Test the object gl is.lengths(gl)
Function to test if an object is of the class multi_tpfit.
is.multi_tpfit(object)is.multi_tpfit(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a multi_tpfit object.
A logical value.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # multidimensional MC models MoPa <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Test the object MoPa is.multi_tpfit(MoPa)data(ACM) # Estimate the parameters of a # multidimensional MC models MoPa <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Test the object MoPa is.multi_tpfit(MoPa)
Function to test if an object is of the class multi_transiogram.
is.multi_transiogram(object)is.multi_transiogram(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a multi_transiogram object.
A logical value.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model TrPr <- predict(RTm, lags) # Test the object TrPr is.multi_transiogram(TrPr)data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model TrPr <- predict(RTm, lags) # Test the object TrPr is.multi_transiogram(TrPr)
The function plots -D sections of a multi-directional transiogram computed without any ellipsoidal interpolation.
is.pemt(object)is.pemt(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a pemt object.
A logical value.
Luca Sartore [email protected]
data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(20, 10, 5), which.dire=c(1, 3), mle = TRUE) # Test the object psEmpTr is.pemt(psEmpTr)data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(20, 10, 5), which.dire=c(1, 3), mle = TRUE) # Test the object psEmpTr is.pemt(psEmpTr)
Function to test if an object is of the class tpfit.
is.tpfit(object)is.tpfit(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a tpfit object.
A logical value.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # one-dimensional MC model MoPa <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Test the object MoPa is.tpfit(MoPa)data(ACM) # Estimate the parameters of a # one-dimensional MC model MoPa <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Test the object MoPa is.tpfit(MoPa)
Function to test if an object is of the class transiogram.
is.transiogram(object)is.transiogram(object)
object |
object to be tested. |
The function returns TRUE if and only if its argument is a transiogram object.
A logical value.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Compute theoretical transition probabilities # from the one-dimensional MC model TTPr <- predict(RTm, lags = 0:2/2) # Compute empirical transition probabilities ETPr <- transiogram(ACM$MAT5, ACM[, 1:3], c(0, 0, 1), 200, 20) # Test the objects TTPr and ETPr is.transiogram(TTPr) is.transiogram(ETPr)data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Compute theoretical transition probabilities # from the one-dimensional MC model TTPr <- predict(RTm, lags = 0:2/2) # Compute empirical transition probabilities ETPr <- transiogram(ACM$MAT5, ACM[, 1:3], c(0, 0, 1), 200, 20) # Test the objects TTPr and ETPr is.transiogram(TTPr) is.transiogram(ETPr)
The function makes a graphical representation of transition probabilities by the use of multiple transiograms.
mixplot(x, main, legend = TRUE, ...)mixplot(x, main, legend = TRUE, ...)
x |
a |
main |
the main title (on top) whose font and size are fixed. |
legend |
a logical value for printing the legend in the graphic. It is |
... |
other arguments to pass to the function |
Transiogram is a diagram which is drawn for a single pair of categories in the direction . It shows the transition probabilities in the -axis for some specific lags in the -axis.
This function permits a graphical approach to compare theoretical vs. empirical transition probabilities for multiple directions.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Li, W. (2007) Transiograms for Characterizing Spatial Variability of Soil Classes. Soil Science Society of America Journal, 71(3), 881-893.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
transiogram, tpfit, predict.tpfit, plot.transiogram, image.multi_tpfit, plot
data(ACM) # Estimate empirical transition # probabilities by points ETr <- transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 100) # Estimate the transition rate matrix RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model TPr <- predict(RTm, lags = ETr$lags) # Plot empirical vs. theoretical transition probabilities mixplot(list(ETr, TPr), type = c("p", "l"), pch = "+", col = c(3, 1))data(ACM) # Estimate empirical transition # probabilities by points ETr <- transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 100) # Estimate the transition rate matrix RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model TPr <- predict(RTm, lags = ETr$lags) # Plot empirical vs. theoretical transition probabilities mixplot(list(ETr, TPr), type = c("p", "l"), pch = "+", col = c(3, 1))
The function estimates the mean length for a -D spatial embedded Markov chain for a specified direction .
mlen(data, coords, loc.id, direction, mle = "avg")mlen(data, coords, loc.id, direction, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
loc.id |
a vector of |
direction |
a |
mle |
a character value. If |
The mean length is the total length occupied by the -th category divided by the number of its embedded occurrences along lines in the direction . More robust methods are implemented, such as the trimmed mean and the trimmed median.
If the stratum lengths are censored, the maximum likelihood approach is more appropriate than the arithmetic mean. In this case, the stratum lengths are assumed to be independent realizations from a log-normal random variable. The quantity to maximize is
where and are vectors of parameters, is the observed stratum length, denotes the upper bound of the censor and denotes a dummy variable which assumes value 1 if and only if the -th stratum is referred to the -th category.
A numeric vector containing the mean length for each observed category.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate the mean lengths for each observed category ml <- mlen(ACM$MAT5, ACM[, 1:3], loc.id, direction, mle = "avg") # Equivalently gl <- getlen(ACM$MAT5, ACM[, 1:3], loc.id, direction, zero.allowed = TRUE) ml1 <- tapply(gl$length, gl$categories, mean)data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate the mean lengths for each observed category ml <- mlen(ACM$MAT5, ACM[, 1:3], loc.id, direction, mle = "avg") # Equivalently gl <- getlen(ACM$MAT5, ACM[, 1:3], loc.id, direction, zero.allowed = TRUE) ml1 <- tapply(gl$length, gl$categories, mean)
The function estimates the model parameters of a -D continuous lag spatial Markov chain. Transition rates matrices along axial directions and proportions of categories are computed.
multi_tpfit(data, coords, method = "ml", tolerance = pi/8, rotation = NULL, max.it = 9000, mle = "avg", ...)multi_tpfit(data, coords, method = "ml", tolerance = pi/8, rotation = NULL, max.it = 9000, mle = "avg", ...)
data |
a categorical data vector of length |
coords |
an |
method |
a character object specifying the method to estimate the transition rates. Possible choises are |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
rotation |
a numerical vector of length |
max.it |
a numerical value which denotes the maximum number of iterations to perform during the optimization phase. It is |
mle |
a character value to pass to the function |
... |
other arguments to pass to the functions |
A -D continuous-lag spatial Markov chain is probabilistic model which is developed by interpolation of the transition rate matrices computed for the main directions. It defines transition probabilities through
where is the lag vector and the entries of are ellipsoidally interpolated.
The ellipsoidal interpolation is given by
where is a standard basis for a -D space.
If the respective entries are replaced by , which is computed as
where and respectively denote the proportions for the -th and -th categories. In so doing, the model may describe the anisotropy of the process.
An object of the class multi_tpfit is returned. The function print.multi_tpfit is used to print the fitted model. The object is a list with the following components:
coordsnames |
a character vector containing the name of each axis. |
coefficients |
a list containing the transition rates matrices computed for each axial direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.multi_tpfit, print.multi_tpfit, image.multi_tpfit, tpfit
data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit(ACM$PERM, ACM[, 1:3])data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit(ACM$PERM, ACM[, 1:3])
The function estimates the model parameters of a -D continuous lag spatial Markov chain by the use of the iterated least squares and the bound-constrained Lagrangian methods. Transition rates matrices along axial directions and proportions of categories are computed.
multi_tpfit_ils(data, coords, max.dist = Inf, mpoints = 20, tolerance = pi/8, rotation = NULL, q = 10, echo = FALSE, ..., mtpfit)multi_tpfit_ils(data, coords, max.dist = Inf, mpoints = 20, tolerance = pi/8, rotation = NULL, q = 10, echo = FALSE, ..., mtpfit)
data |
a categorical data vector of length |
coords |
an |
max.dist |
a numerical value which defines the maximum lag value. It is |
mpoints |
a numerical value which defines the number of lag intervals. |
tolerance |
a numerical value for the tolerance angle (in radians). It is |
rotation |
a numerical vector of length |
q |
a numerical value greater than one for a constant which controls the growth of the penalization term in the loss function. It is equal to |
echo |
a logical value; if |
... |
other arguments to pass to the function |
mtpfit |
an object |
A -D continuous-lag spatial Markov chain is probabilistic model which is developed by interpolation of the transition rate matrices computed for the main directions. It defines transition probabilities through
where is the lag vector and the entries of are ellipsoidally interpolated.
The ellipsoidal interpolation is given by
where is a standard basis for a -D space.
If the respective entries are replaced by , which is computed as
where and respectively denote the proportions for the -th and -th categories. In so doing, the model may describe the anisotropy of the process.
In particular, to estimate entries of transition rate matrices computed for the main axial directions, we need to minimize the discrepancies between the empirical transiograms (see transiogram) and the predicted transition probabilities.
By the use of the iterated least squares, the diagonal entries of are constrained to be negative, while the off-diagonal transition rates are constrained to be positive. Further constraints are considered in order to obtain a proper transition rates matrix.
An object of the class multi_tpfit is returned. The function print.multi_tpfit is used to print the fitted model. The object is a list with the following components:
coordsnames |
a character vector containing the name of each axis. |
coefficients |
a list containing the transition rates matrices computed for each axial direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
If the process is not stationary, the optimization algorithm does not converge.
Luca Sartore [email protected]
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.multi_tpfit, print.multi_tpfit, image.multi_tpfit, tpfit_ils, transiogram
data(ACM) # Estimate the parameters of a # multidimensional MC model multi_tpfit_ils(ACM$MAT3, ACM[, 1:3], 100)data(ACM) # Estimate the parameters of a # multidimensional MC model multi_tpfit_ils(ACM$MAT3, ACM[, 1:3], 100)
The function estimates the model parameters of a -D continuous lag spatial Markov chain. Transition rates matrices along axial directions and proportions of categories are computed.
multi_tpfit_me(data, coords, tolerance = pi/8, max.it = 9000, rotation = NULL, mle = "avg")multi_tpfit_me(data, coords, tolerance = pi/8, max.it = 9000, rotation = NULL, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
tolerance |
a numerical value for the tolerance angle (in radians). It is |
max.it |
a numerical value which denotes the maximum number of iterations to perform during the optimization phase. It is |
rotation |
a numerical vector of length |
mle |
a character value to pass to the function |
A -D continuous-lag spatial Markov chain is probabilistic model which is developed by interpolation of the transition rate matrices computed for the main directions by the use of the function tpfit_me. It defines transition probabilities through
where is the lag vector and the entries of are ellipsoidally interpolated.
The ellipsoidal interpolation is given by
where is a standard basis for a -D space.
If the respective entries are replaced by , which is computed as
where and respectively denote the proportions for the -th and -th categories. In so doing, the model may describe the anisotropy of the process.
When some entries of the rates matrices are not identifiable, it is suggested to vary the tolerance coefficient and the rotation angles. This problem may be also avoided if the input argument mle is to set to be "mlk".
An object of the class multi_tpfit is returned. The function print.multi_tpfit is used to print the fitted model. The object is a list with the following components:
coordsnames |
a character vector containing the name of each axis. |
coefficients |
a list containing the transition rates matrices computed for each axial direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
predict.multi_tpfit, print.multi_tpfit, image.multi_tpfit, tpfit_me
data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit_me(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit_me(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit_me(ACM$PERM, ACM[, 1:3])data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit_me(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit_me(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit_me(ACM$PERM, ACM[, 1:3])
The function estimates the model parameters of a -D continuous lag spatial Markov chain. Transition rates matrices along axial directions and proportions of categories are computed.
multi_tpfit_ml(data, coords, tolerance = pi/8, rotation = NULL, mle = "avg")multi_tpfit_ml(data, coords, tolerance = pi/8, rotation = NULL, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
rotation |
a numerical vector of length |
mle |
a character value to pass to the function |
A -D continuous-lag spatial Markov chain is probabilistic model which is developed by interpolation of the transition rate matrices computed for the main directions. It defines transition probabilities through
where is the lag vector and the entries of are ellipsoidally interpolated.
The ellipsoidal interpolation is given by
where is a standard basis for a -D space.
If the respective entries are replaced by , which is computed as
where and respectively denote the proportions for the -th and -th categories. In so doing, the model may describe the anisotropy of the process.
When some entries of the rates matrices are not identifiable, it is suggested to vary the tolerance coefficient and the rotation angles. This problem may be also avoided if the input argument mle is to set to be "mlk".
An object of the class multi_tpfit is returned. The function print.multi_tpfit is used to print the fitted model. The object is a list with the following components:
coordsnames |
a character vector containing the name of each axis. |
coefficients |
a list containing the transition rates matrices computed for each axial direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.multi_tpfit, print.multi_tpfit, image.multi_tpfit, tpfit_ml
data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit_ml(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit_ml(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit_ml(ACM$PERM, ACM[, 1:3])data(ACM) # Estimate transition rates matrices and # proportions for the categorical variable MAT5 multi_tpfit_ml(ACM$MAT5, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable MAT3 multi_tpfit_ml(ACM$MAT3, ACM[, 1:3]) # Estimate transition rates matrices and # proportions for the categorical variable PERM multi_tpfit_ml(ACM$PERM, ACM[, 1:3])
The function computes the multi-directional transiograms without any ellipsoidal interpolation for -D sections.
pemt(data, coords, mpoints, which.dire, max.dist, tolerance = pi/8, rotation = NULL, mle = "avg")pemt(data, coords, mpoints, which.dire, max.dist, tolerance = pi/8, rotation = NULL, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
mpoints |
the number of points per axes. It controls the accuracy of images to plot. |
which.dire |
a vector with two chosen axial directions. If omitted, all |
max.dist |
a scalar or a vector of maximum length for the chosen axial directions. |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
rotation |
a numerical vector of length |
mle |
a character value to pass to the function |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. The probability is computed for any lag vector through
where entries of are not ellipsoidally interpolated, but they are estimated for the direction specified by the vector .
In particular cases, some entries of the estimated matrix might be not finite, so that the exponential matrix is computable and the resulting transition probabilities are set to be NaN. If mle = "mlk", this problem may be partially solved.
The exponential matrix is evaluated by the scaling and squaring algorithm.
An object of class pemt is returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
multi_tpfit_ml, tpfit_ml, image.pemt, plot.transiogram
data(ACM) # Compute a 2-D section of a # multi-directional transiogram pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "mdn")data(ACM) # Compute a 2-D section of a # multi-directional transiogram pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire=c(1, 3), mle = "mdn")
The function draws perspective-plots the -D sections of a predicted multidimensional transiograms computed through ellipsoidal interpolation.
## S3 method for class 'multi_tpfit' persp(x, mpoints, which.dire, max.dist, main, mar, ask = TRUE, col = "white", ...)## S3 method for class 'multi_tpfit' persp(x, mpoints, which.dire, max.dist, main, mar, ask = TRUE, col = "white", ...)
x |
an object of the class |
mpoints |
the number of points per axes. It controls the accuracy of images to plot. |
which.dire |
a vector with two chosen axial directions. If omitted, all |
max.dist |
a scalar or a vector of maximum length for the chosen axial directions. |
main |
the main title (on top) whose font and size are fixed. |
mar |
a scalar or a numerical vector of the form |
ask |
a logical value; if |
col |
a list of colors which is usually generated by |
... |
other arguments to pass to the function |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. It is computed for any lag vector through
where entries of are ellipsoidally interpolated (see multi_tpfit).
The exponential matrix is evaluated by the scaling and squaring algorithm.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
multi_tpfit, persp.multi_tpfit, persp, pemt, persp.pemt, plot.transiogram
data(ACM) # Estimate model parameter x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Set short names for categories 3 and 4 names(x$prop)[3:4] <- c("Clay and Sand", "Gravel and Sand") # 3D-Plot for a 2-D theoretical sections of # a multidimensional transiogram persp(x, 15, max.dist = c(200, 200, 20), which.dire = 2:3, mar = .7, col = rainbow(500), theta = 15, phi = 45)data(ACM) # Estimate model parameter x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Set short names for categories 3 and 4 names(x$prop)[3:4] <- c("Clay and Sand", "Gravel and Sand") # 3D-Plot for a 2-D theoretical sections of # a multidimensional transiogram persp(x, 15, max.dist = c(200, 200, 20), which.dire = 2:3, mar = .7, col = rainbow(500), theta = 15, phi = 45)
The function draws perspective-plots the -D sections of a multi-directional transiogram computed without any ellipsoidal interpolation.
## S3 method for class 'pemt' persp(x, main, mar, ask = TRUE, col = "white", ...)## S3 method for class 'pemt' persp(x, main, mar, ask = TRUE, col = "white", ...)
x |
an object of the class |
main |
the main title (on top) whose font and size are fixed. |
mar |
a scalar or a numerical vector of the form |
ask |
a logical value; if |
col |
a list of colors which is usually generated by |
... |
other arguments to pass to the function |
A multidimensional transiogram is a diagram which shows the transition probabilities for a single pair of categories. The probability is computed for any lag vector through
where entries of are not ellipsoidally interpolated, but they are estimated for the direction specified by the vector .
The exponential matrix is evaluated by the scaling and squaring algorithm.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Higham, N. J. (2008) Functions of Matrices: Theory and Computation. Society for Industrial and Applied Mathematics.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
pemt, persp.multi_tpfit, persp, multi_tpfit, image.pemt, plot.transiogram
data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire = c(1, 3)) # 3D-Plot for a 2-D sections of # a multi-directional transiogram persp(psEmpTr, col = rainbow(500), mar = .7, theta = 15, phi = 45)data(ACM) # Compute a 2-D section of a # multi-directional transiogram psEmpTr <- pemt(ACM$MAT3, ACM[, 1:3], 2, max.dist = c(200, 200, 20), which.dire = c(1, 3)) # 3D-Plot for a 2-D sections of # a multi-directional transiogram persp(psEmpTr, col = rainbow(500), mar = .7, theta = 15, phi = 45)
The function plot the empirical densities of stratum lengths computed along a given direction.
## S3 method for class 'density.lengths' plot(x, main = NULL, xlab = NULL, ylab = "Density", type = "l", zero.line = TRUE, ...)## S3 method for class 'density.lengths' plot(x, main = NULL, xlab = NULL, ylab = "Density", type = "l", zero.line = TRUE, ...)
x |
an object of the class |
main |
an overall title for the plot. |
xlab |
a title for the |
ylab |
a title for the |
type |
plotting parameter for the type of graphic (see |
zero.line |
logical value. If |
... |
other plotting parameters. |
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
density.default, density.lengths, plot, print.density.lengths
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum log-lengths dgl <- density(gl, log = TRUE) # Plot the empirical densities of stratum log-lengths plot(dgl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum log-lengths dgl <- density(gl, log = TRUE) # Plot the empirical densities of stratum log-lengths plot(dgl)
The function plots objects of class hist.lengths.
## S3 method for class 'hist.lengths' plot(x, ...)## S3 method for class 'hist.lengths' plot(x, ...)
x |
an object of the class |
... |
further plotting parameters. |
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
hist, hist.lengths, plot, print.density.lengths
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the histograms hgl <- hist(gl, plot = FALSE) # Plot the histograms plot(hgl, col = "#efffef")data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the histograms hgl <- hist(gl, plot = FALSE) # Plot the histograms plot(hgl, col = "#efffef")
The function makes a graphical representation of the stratum lengths.
## S3 method for class 'lengths' plot(x, ..., log = FALSE, zeros.rm = TRUE)## S3 method for class 'lengths' plot(x, ..., log = FALSE, zeros.rm = TRUE)
x |
an object of the class |
... |
other arguments to pass to the function |
log |
a logical value. If |
zeros.rm |
a logical value. If |
The box-and-whisker plots give some information about the distribution of the stratum lengths for the observed categories along a given direction.
An image is produced on the current graphics device; by the use of boxplot.lengths, the same image is produced. The function returns a list with the following components:
stats |
a matrix containing the values used to plot the box-and-whisker plots. |
n |
a vector with the number of observations for each category. |
conf |
a matrix containing further values to draw the lower and upper extremes of the notch. |
out |
a vectors with the values of the outlier points. |
group |
a vector whose elements indicate to which category the outlier belongs. |
names |
a character vector with the names of each category. |
Luca Sartore [email protected]
boxplot.lengths, boxplot, getlen
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Plot the object gl plot(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Plot the object gl plot(gl)
The function makes a graphical representation of transition probabilities by the use of transiogram.
## S3 method for class 'transiogram' plot(x, ..., main, legend = FALSE, ci = NULL)## S3 method for class 'transiogram' plot(x, ..., main, legend = FALSE, ci = NULL)
x |
an object of the class |
... |
other arguments to pass to the function |
main |
the main title (on top) whose font and size are fixed. |
legend |
a logical value; if |
ci |
a numerical value in the interval (0, 1) denoting the confidence of the interval around transition probabilities. If |
Transiogram is a diagram which is drawn for a single pair of categories in the direction . It shows the transition probabilities in the -axis for some specific lags in the -axis.
Confidence intervals are computed on the log odds of the transition probabilities. The approximation of the confidence bounds is based on the delta method applied on the logistic transformation.
An image is produced on the current graphics device. No values are returned.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Li, W. (2007) Transiograms for Characterizing Spatial Variability of Soil Classes. Soil Science Society of America Journal, 71(3), 881-893.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
tpfit, predict.tpfit, mixplot, image.multi_tpfit, plot
data(ACM) # Estimate empirical transition # probabilities by points ETr <- transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 100, 100) # Estimate the transition rate matrix RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model TPr <- predict(RTm, lags = ETr$lags) # Plot empirical transition probabilities plot(ETr, type = "l", ci = 0.99) # Plot theoretical transition probabilities plot(TPr, type = "l")data(ACM) # Estimate empirical transition # probabilities by points ETr <- transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 100, 100) # Estimate the transition rate matrix RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model TPr <- predict(RTm, lags = ETr$lags) # Plot empirical transition probabilities plot(ETr, type = "l", ci = 0.99) # Plot theoretical transition probabilities plot(TPr, type = "l")
The function computes theoretical transition probabilities of a -D continuous-lag spatial Markov chain for a specified set of lags.
## S3 method for class 'multi_tpfit' predict(object, lags, byrow = TRUE, ...)## S3 method for class 'multi_tpfit' predict(object, lags, byrow = TRUE, ...)
object |
an object of the class |
lags |
a lag vector or matrix of |
byrow |
a logical value; if |
... |
further arguments passed from other methods. |
A -D continuous-lag spatial Markov chain is probabilistic model which is developed by interpolation of the transition rate matrices computed for the main directions. It defines the transition probability through the entry of the following matrix
where is the lag vector and the entries of are ellipsoidally interpolated.
An object of the class multi_transiogram is returned. The print.multi_transiogram function is used to print computed probabilities. The object is a list with the following components:
Tmat |
a 3-D array containing the probabilities. |
lags |
a matrix containing the lag vectors. |
type |
a character string which specifies that computed probabilities are theoretical. |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
multi_tpfit, print.multi_tpfit, image.multi_tpfit, tpfit, transiogram
data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model predict(RTm, lags)data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model predict(RTm, lags)
The function computes theoretical transition probabilities of a 1-D continuous-lag spatial Markov chain for a specified set of lags.
## S3 method for class 'tpfit' predict(object, lags, ...)## S3 method for class 'tpfit' predict(object, lags, ...)
object |
an object of the class |
lags |
a vector of 1-D lags. |
... |
further arguments passed from other methods. |
A 1-D continuous-lag spatial Markov chain is probabilistic model which involves a transition rate matrix computed for the direction . It defines the transition probability through the entry of the following matrix
where is a positive lag value.
An object of the class transiogram is returned. The function print.transiogram is used to print computed probabilities. The object is a list with the following components:
Tmat |
a 3-D array containing the probabilities. |
lags |
a vector containing one-dimensional lags. |
type |
a character string which specifies that computed probabilities are theoretical. |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
tpfit, print.tpfit, plot.transiogram, transiogram, multi_tpfit
data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model predict(RTm, lags = 0:2/2)data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT3, ACM[, 1:3], c(0, 0, 1)) # Compute transition probabilities # from the one-dimensional MC model predict(RTm, lags = 0:2/2)
he function a summary of the empirical density stratum lengths calculated by density.lengths.
## S3 method for class 'density.lengths' print(x, digits = NULL, ...)## S3 method for class 'density.lengths' print(x, digits = NULL, ...)
x |
an object of the class |
digits |
minimal number of digits, see |
... |
further arguments to pass to the function |
A summary of the empirical distributions is printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
density.lengths, plot.density.lengths
data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum lengths dgl <- density(gl) # Print the empirical densities of stratum lengths print(dgl)data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Compute the empirical densities of stratum lengths dgl <- density(gl) # Print the empirical densities of stratum lengths print(dgl)
The function prints stratum lengths given by getlen.
## S3 method for class 'lengths' print(x, ...)## S3 method for class 'lengths' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
Stratum lengths grouped by category are printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Print stratum lengths print(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Print stratum lengths print(gl)
The function prints parameter estimation results given by multi_tpfit.
## S3 method for class 'multi_tpfit' print(x, ...)## S3 method for class 'multi_tpfit' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
Estimation results are printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # multidimensional MC models MoPa <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Print results print(MoPa)data(ACM) # Estimate the parameters of a # multidimensional MC models MoPa <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Print results print(MoPa)
The function prints theoretical transition probabilities given by predict.multi_tpfit.
## S3 method for class 'multi_transiogram' print(x, ...)## S3 method for class 'multi_transiogram' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
Transition probabilities are printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model TrPr <- predict(RTm, lags) # Print results print(TrPr)data(ACM) # Estimate the parameters of a # multidimensional MC model RTm <- multi_tpfit(ACM$MAT3, ACM[, 1:3]) # Generate the matrix of # multidimensional lags lags <- expand.grid(X=-1:1, Y=-1:1, Z=-1:1) lags <- as.matrix(lags) # Compute transition probabilities # from the multidimensional MC model TrPr <- predict(RTm, lags) # Print results print(TrPr)
The function prints the summary of stratum lengths given by summary.lengths.
## S3 method for class 'summary.lengths' print(x, ...)## S3 method for class 'summary.lengths' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
The summary of stratum lengths grouped by category is printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Summarize the stratum lengths sgl <- summary(gl) # Print the summary of stratum lengths print(sgl)data(ACM) direction <- c(0,0,1) # Compute the appartaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Summarize the stratum lengths sgl <- summary(gl) # Print the summary of stratum lengths print(sgl)
The function prints parameter estimation results given by tpfit.
## S3 method for class 'tpfit' print(x, ...)## S3 method for class 'tpfit' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
Estimation results are printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # one-dimensional MC model MoPa <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Print results print(MoPa)data(ACM) # Estimate the parameters of a # one-dimensional MC model MoPa <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Print results print(MoPa)
The function prints transition probabilities given by predict.multi_tpfit or transiogram.
## S3 method for class 'transiogram' print(x, ...)## S3 method for class 'transiogram' print(x, ...)
x |
an object of the class |
... |
further arguments passed to or from other methods. |
Transition probabilities are printed on the screen or other output devices. No values are returned.
Luca Sartore [email protected]
data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Compute theoretical transition probabilities # from the one-dimensional MC model TTPr <- predict(RTm, lags = 0:2/2) # Compute empirical transition probabilities ETPr <- transiogram(ACM$MAT5, ACM[, 1:3], c(0, 0, 1), 200, 20) # Print results print(TTPr) print(ETPr)data(ACM) # Estimate the parameters of a # one-dimensional MC model RTm <- tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1)) # Compute theoretical transition probabilities # from the one-dimensional MC model TTPr <- predict(RTm, lags = 0:2/2) # Compute empirical transition probabilities ETPr <- transiogram(ACM$MAT5, ACM[, 1:3], c(0, 0, 1), 200, 20) # Print results print(TTPr) print(ETPr)
The function adjusts a simulated random field generated by the sim function.
quench(x, data, coords, sim, GA = FALSE, optype = c("param", "fullprobs", "semiprobs", "coordprobs"), max.it = 1000, knn = 12)quench(x, data, coords, sim, GA = FALSE, optype = c("param", "fullprobs", "semiprobs", "coordprobs"), max.it = 1000, knn = 12)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
sim |
an object of the class |
GA |
a logical value; if |
optype |
a character which denotes the objective function to compute when the optimization is performed. |
max.it |
a numerical value which specifies the maximum number of iterations to stop the optimization algorithm. For proper results, it should be a multiple of the number of simulation points. |
knn |
an integer value which specifies the number of k-nearest neighbours for each simulation point. An optimal number is between 4 and 12. If |
This method perform a simulated annealing or a genetic algorithm to modify the simulation results, in order to reduce artifacts effects. In practice, each simulated configuration is adjusted to reach a pattern similar to the observed sample data. There are several objective functions for this purpose, by setting optype equal to "param" the optimization is performed through parametric methods. The alternatives "fullprobs" and "semiprobs" are based on transition probabilities computed among simulation points, while the option "coordprobs" is based on transition probabilities calculated among observation and simulation points.
This procedure should be executed by setting max.it equal at least to the simulation grid size, or its multiples.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1996) Transition Probability-Based Indicator Geostatistics. Mathematical Geosciences, 28(4), 453-476.
Carle, S. F. (1999) T-PROGS: Transition Probability Geostatistical Software. University of California, Davis.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
Weise, T. (2009) Global Optimization Algorithms - Theory and Application. https://archive.org/details/Thomas_Weise__Global_Optimization_Algorithms_Theory_and_Application.
sim_ck, sim_ik, sim_mcs, sim_path
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Ordinary Indicator Kriging algorithm myOIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid) # Perform the quenching algorithm # to adjust simulation quench(x, ACM$MAT5, ACM[, 1:3], myOIKSim, optype = "coordprobs", max.it = 2, knn = 12)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Ordinary Indicator Kriging algorithm myOIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid) # Perform the quenching algorithm # to adjust simulation quench(x, ACM$MAT5, ACM[, 1:3], myOIKSim, optype = "coordprobs", max.it = 2, knn = 12)
The function set the number of CPU cores for parallel computation by the use of OpenMP library (https://www.openmp.org/). If the package was not complied with the library OpenMP (>= 3.0), this function is disabled.
setCores(n)setCores(n)
n |
an integer value denoting the number of CPU cores to use; if it exceeds the total number of cores, all of them will be used. If missing, the number of CPU cores in use will be displayed. |
When the package is loaded, only one CPU core is used.
The total number of CPU cores in use will be returned and a message will be displayed. If the package was not complied with the library OpenMP (>= 3.0), the value one will be returned.
Luca Sartore [email protected]
SunTM ONE Studio 8 (2003) OpenMP API User's Guide. Sun Microsystems Inc., Santa Clara, U.S.A.
#Display the number of CPU cores in use setCores() #Set 2 CPU cores for parallel computation setCores(2) #Set 1 CPU core for serial computation setCores(1)#Display the number of CPU cores in use setCores() #Set 2 CPU cores for parallel computation setCores(2) #Set 1 CPU core for serial computation setCores(1)
The function simulates a random field. The simulation methods available are based on Indicator Kriging techniques (IK and CK), Fixed and Random Path (PATH) and Multinomial Categorical Simulation (MCS).
sim(x, data, coords, grid, method = "ik", ..., entropy = FALSE)sim(x, data, coords, grid, method = "ik", ..., entropy = FALSE)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
grid |
an |
method |
a character object specifying the method to simulate the random field. Possible choises are |
... |
other arguments to pass to the functions |
entropy |
a logical value. If |
The methods implemented compute the approximation of posterior probabilities
Once the probabilities are calculated for all the points in the simulation grid, the predictions (based on most probable category) and simulations are returned.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities is returned. Two extra columns respectively denoting the entropy and standardized entorpy are bindend to the data frame when argument entropy = TRUE.
Allard, D., D'Or, D., Froidevaux, R. (2011) An efficient maximum entropy approach for categorical variable prediction. European Journal of Soil Science, 62(3), 381-393.
Carle, S. F., Fogg, G. E. (1996) Transition Probability-Based Indicator Geostatistics. Mathematical Geosciences, 28(4), 453-476.
Carle, S. F. (1999) T-PROGS: Transition Probability Geostatistical Software. University of California, Davis.
Li, W. (2007) A Fixed-Path Markov Chain Algorithm for Conditional Simulation of Discrete Spatial Variables. Mathematical Geology, 39(2), 159-176.
Li, W. (2007) Markov Chain Random Fields for Estimation of Categorical Variables. Mathematical Geology, 39(June), 321-335.
Pickard, D. K. (1980) Unilateral Markov Fields. Advances in Applied Probability, 12(3), 655-671.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
Weise, T. (2009) Global Optimization Algorithms - Theory and Application. https://archive.org/details/Thomas_Weise__Global_Optimization_Algorithms_Theory_and_Application.
sim_ik, sim_ck, sim_path, sim_mcs
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Kriging algorithm and mySim <- sim(x, ACM$MAT5, ACM[, 1:3], mygrid)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Kriging algorithm and mySim <- sim(x, ACM$MAT5, ACM[, 1:3], mygrid)
The function simulates a random field through the Indicator Cokriging technique.
sim_ck(x, data, coords, grid, knn = 12, ordinary = TRUE, entropy = FALSE)sim_ck(x, data, coords, grid, knn = 12, ordinary = TRUE, entropy = FALSE)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
grid |
an |
knn |
an integer value which specifies the number of k-nearest neighbours for each simulation point. An optimal number is between 4 and 12. If |
ordinary |
a logical value; if |
entropy |
a logical value. If |
This method computes an approximation of posterior probabilities
The probability is calculated as the weighted sum of indicator variables which denote the presence of the -th category in observed points . Weights involved in the sum are the solution of a system of equations.
Probabilities approximated are usually truncated and normalized with respect to the probability constraints, because such probabilities might lie outside the interval . The normalization procedure is designed such that it is not possible to obtain vectors such that the sum of their probabilities is always equal to one.
When an initial configuration is simulated, it should be modified to reach a pattern similar to the sample by the use of the quench function.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities is returned. Two extra columns respectively denoting the entropy and standardized entorpy are bindend to the data frame when argument entropy = TRUE.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1996) Transition Probability-Based Indicator Geostatistics. Mathematical Geosciences, 28(4), 453-476.
Carle, S. F. (1999) T-PROGS: Transition Probability Geostatistical Software. University of California, Davis.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
Weise, T. (2009) Global Optimization Algorithms - Theory and Application. https://archive.org/details/Thomas_Weise__Global_Optimization_Algorithms_Theory_and_Application.
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Cokriging algorithm mySCKSim <- sim_ck(x, ACM$MAT5, ACM[, 1:3], mygrid, ordinary = FALSE) # Simulate the random field through # Ordinary Indicator Cokriging algorithm myOCKSim <- sim_ck(x, ACM$MAT5, ACM[, 1:3], mygrid)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Cokriging algorithm mySCKSim <- sim_ck(x, ACM$MAT5, ACM[, 1:3], mygrid, ordinary = FALSE) # Simulate the random field through # Ordinary Indicator Cokriging algorithm myOCKSim <- sim_ck(x, ACM$MAT5, ACM[, 1:3], mygrid)
The function simulates a random field through the Indicator Kriging technique.
sim_ik(x, data, coords, grid, knn = 12, ordinary = TRUE, entropy = FALSE)sim_ik(x, data, coords, grid, knn = 12, ordinary = TRUE, entropy = FALSE)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
grid |
an |
knn |
an integer value which specifies the number of k-nearest neighbours for each simulation point. An optimal number is between 4 and 12. If |
ordinary |
a logical value; if |
entropy |
a logical value. If |
This method computes an approximation of posterior probabilities
The probability is calculated as the sum of the observed proportion and the weighted sum of indicator variables which denote the presence of the -th category in observed points . Weights involved in the sum are the solution of a system of equations.
Probabilities approximated are usually truncated and normalized with respect to the probability constraints, because such probabilities might lie outside the interval . The normalization procedure is designed such that it is not possible to obtain vectors such that the sum of their probabilities is always equal to one.
When an initial configuration is simulated, it should be modified to reach a pattern similar to the sample by the use of the quench function.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities is returned. Two extra columns respectively denoting the entropy and standardized entorpy are bindend to the data frame when argument entropy = TRUE.
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1996) Transition Probability-Based Indicator Geostatistics. Mathematical Geosciences, 28(4), 453-476.
Carle, S. F. (1999) T-PROGS: Transition Probability Geostatistical Software. University of California, Davis.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
Weise, T. (2009) Global Optimization Algorithms - Theory and Application. https://archive.org/details/Thomas_Weise__Global_Optimization_Algorithms_Theory_and_Application.
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Kriging algorithm mySIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid, ordinary = FALSE) # Simulate the random field through # Ordinary Indicator Kriging algorithm myOIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # Simple Indicator Kriging algorithm mySIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid, ordinary = FALSE) # Simulate the random field through # Ordinary Indicator Kriging algorithm myOIKSim <- sim_ik(x, ACM$MAT5, ACM[, 1:3], mygrid)
The function simulates a random field through the Multinomial Categorical Simulation technique (MCS).
sim_mcs(x, data, coords, grid, knn = NULL, entropy = FALSE)sim_mcs(x, data, coords, grid, knn = NULL, entropy = FALSE)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
grid |
an |
knn |
an integer value which specifies the number of k-nearest neighbours for each simulation point. If |
entropy |
a logical value. If |
This method computes an approximation of posterior probabilities
The algorithm is based on the Bayesian maximum entropy approach and it honours both the model structure and observed data.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities is returned. Two extra columns respectively denoting the entropy and standardized entorpy are bindend to the data frame when argument entropy = TRUE.
Luca Sartore [email protected]
Allard, D., D'Or, D., Froidevaux, R. (2011) An efficient maximum entropy approach for categorical variable prediction. European Journal of Soil Science, 62(3), 381-393.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 3) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 3) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field myMCSim <- sim_mcs(x, ACM$MAT5, ACM[, 1:3], mygrid)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 3) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 3) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field myMCSim <- sim_mcs(x, ACM$MAT5, ACM[, 1:3], mygrid)
The function simulates a random field through the Fixed Path algorithm or Random Path technique.
sim_path(x, data, coords, grid, radius, fixed = FALSE, entropy = FALSE)sim_path(x, data, coords, grid, radius, fixed = FALSE, entropy = FALSE)
x |
an object of the class |
data |
a categorical data vector of length |
coords |
an |
grid |
an |
radius |
a numerical value that specifies a proper radius to search the nearest observed points within a |
fixed |
a logical value; if |
entropy |
a logical value. If |
These methods compute an approximation of posterior probabilities
Path algorithms are based on Pickard random fields, so that the states of such chain at any unsampled location depends on the state of its nearest known neighbours in axial directions.
A data frame containing the simulation grid, the simulated random field, predicted values and the approximated probabilities is returned. Two extra columns respectively denoting the entropy and standardized entorpy are bindend to the data frame when argument entropy = TRUE.
Luca Sartore [email protected]
Li, W. (2007) A Fixed-Path Markov Chain Algorithm for Conditional Simulation of Discrete Spatial Variables. Mathematical Geology, 39(2), 159-176.
Li, W. (2007) Markov Chain Random Fields for Estimation of Categorical Variables. Mathematical Geology, 39(June), 321-335.
Pickard, D. K. (1980) Unilateral Markov Fields. Advances in Applied Probability, 12(3), 655-671.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # the fixed path algorithm myFixPathSim <- sim_path(x, ACM$MAT5, ACM[, 1:3], mygrid, radius = 50, fixed = TRUE) # Simulate the random field through # the random path algorithm myRndPathSim <- sim_path(x, ACM$MAT5, ACM[, 1:3], mygrid, radius = 50)data(ACM) # Model parameters estimation for the # multinomial categorical simulation x <- multi_tpfit(ACM$MAT5, ACM[, 1:3]) # Generate the simulation grid mygrid <- list() mygrid$X <- seq(min(ACM$X), max(ACM$X), length = 20) mygrid$Y <- seq(min(ACM$Y), max(ACM$Y), length = 20) mygrid$Z <- -40 * 0:9 - 1 mygrid <- as.matrix(expand.grid(mygrid$X, mygrid$Y, mygrid$Z)) # Simulate the random field through # the fixed path algorithm myFixPathSim <- sim_path(x, ACM$MAT5, ACM[, 1:3], mygrid, radius = 50, fixed = TRUE) # Simulate the random field through # the random path algorithm myRndPathSim <- sim_path(x, ACM$MAT5, ACM[, 1:3], mygrid, radius = 50)
The function summarizes the stratum lengths for each observed category.
## S3 method for class 'lengths' summary(object, ..., zeros.rm = TRUE)## S3 method for class 'lengths' summary(object, ..., zeros.rm = TRUE)
object |
an object of the class |
... |
further arguments passed to or from other methods. |
zeros.rm |
a logical values. If |
An object of class summary.lengths containing the minimum, the first quartile, the median, the mean, the third quartile and the maximum of the stratum lengths for each observed category.
Luca Sartore [email protected]
data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Summarize the stratum lengths sgl <- summary(gl)data(ACM) direction <- c(0,0,1) # Compute the appertaining directional line for each location loc.id <- which_lines(ACM[, 1:3], direction) # Estimate stratum lengths gl <- getlen(ACM$MAT3, ACM[, 1:3], loc.id, direction) # Summarize the stratum lengths sgl <- summary(gl)
The function estimates the model parameters of a 1-D continuous lag spatial Markov chain. Transition rates matrix along a user defined direction and proportions of categories are computed.
tpfit(data, coords, direction, method = "ml", tolerance = pi/8, max.it = 9000, mle = "avg", ...)tpfit(data, coords, direction, method = "ml", tolerance = pi/8, max.it = 9000, mle = "avg", ...)
data |
a categorical data vector of length |
coords |
an |
direction |
a |
method |
a character object specifying the method to estimate the transition rates. Possible choises are |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
max.it |
a numerical value which denotes the maximum number of iterations to perform during the optimization phase. It is |
mle |
a character value to pass to the function |
... |
other arguments to pass to the functions |
A 1-D continuous-lag spatial Markov chain is probabilistic model which involves a transition rate matrix computed for the direction . It defines the transition probability through the entry of the following matrix
where is a positive lag value.
Three methods are available to calculate entries of the transition rate matrix. The mean length method is performed by the use of the function tpfit_ml, the iterated least squares are applied through the function tpfit_ils, while the function tpfit_me implements the maximum entropy method.
An object of the class tpfit is returned. The function print.tpfit is used to print the fitted model. The object is a list with the following components:
coefficients |
the transition rates matrix computed for the user defined direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.tpfit, print.tpfit, multi_tpfit, transiogram
data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1))data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit(ACM$MAT5, ACM[, 1:3], c(0, 0, 1))
The function estimates the model parameters of a 1-D continuous lag spatial Markov chain by the use of the iterated least squares and the bound-constrained Lagrangian methods. Transition rates matrix along a user defined direction and proportions of categories are computed.
tpfit_ils(data, coords, direction, max.dist = Inf, mpoints = 20, tolerance = pi/8, q = 10, echo = FALSE, ..., tpfit)tpfit_ils(data, coords, direction, max.dist = Inf, mpoints = 20, tolerance = pi/8, q = 10, echo = FALSE, ..., tpfit)
data |
a categorical data vector of length |
coords |
an |
direction |
a |
max.dist |
a numerical value which defines the maximum lag value. It's |
mpoints |
a numerical value which defines the number of lag intervals. |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
q |
a numerical value greater than one for a constant which controls the growth of the penalization term in the loss function. It is equal to |
echo |
a logical value; if |
... |
other arguments to pass to the function |
tpfit |
an object |
A 1-D continuous-lag spatial Markov chain is probabilistic model which involves a transition rate matrix computed for the direction . It defines the transition probability through the entry of the following matrix
where is a positive lag value.
To calculate entries of the transition rate matrix, we need to minimize the discrepancies between the empirical transiogram (see transiogram) and the predicted transition probabilities.
By the use of the iterated least squares, the diagonal entries of are constrained to be negative,
while the off-diagonal transition rates are constrained to be positive. Further constraints are considered in order to obtain a proper transition rates matrix.
An object of the class tpfit is returned. The function print.tpfit is used to print the fitted model. The object is a list with the following components:
coefficients |
the transition rates matrix computed for the user defined direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
If the process is not stationary, the optimization algorithm does not converge.
Luca Sartore [email protected]
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.tpfit, print.tpfit, multi_tpfit_ils, transiogram
data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_ils(ACM$MAT3, ACM[, 1:3], c(0,0,1), 100)data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_ils(ACM$MAT3, ACM[, 1:3], c(0,0,1), 100)
The function estimates the model parameters of a 1-D continuous lag spatial Markov chain by the use of the maximum entropy method. Transition rates matrix along a user defined direction and proportions of categories are computed.
tpfit_me(data, coords, direction, tolerance = pi/8, max.it = 9000, mle = "avg")tpfit_me(data, coords, direction, tolerance = pi/8, max.it = 9000, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
direction |
a |
tolerance |
a numerical value for the tolerance angle (in radians). It is |
max.it |
a numerical value which denotes the maximum number of iterations to perform during the optimization phase. It is |
mle |
a character value to pass to the function |
A 1-D continuous-lag spatial Markov chain is probabilistic model which involves a transition rate matrix computed for the direction . It defines the transition probability through the entry of the following matrix
where is a positive lag value.
To calculate entries of the transition rate matrix, we need to maximize the entropy of the transition probabilities of embedded occurrences along a given direction . The entropy is defined as
where are transition probabilities of embedded occurrences. It is maximized by the use of the iterative proportion fitting method.
When some entries of the matrix are not identifiable, it is suggested to vary the tolerance coefficient or to set the input argument mle to "mlk".
An object of the class tpfit is returned. The function print.tpfit is used to print the fitted model. The object is a list with the following components:
coefficients |
the transition rates matrix computed for the user defined direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
predict.tpfit, print.tpfit, multi_tpfit_me
data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_me(ACM$MAT5, ACM[, 1:3], c(0,0,1))data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_me(ACM$MAT5, ACM[, 1:3], c(0,0,1))
The function estimates the model parameters of a 1-D continuous lag spatial Markov chain by the use of the mean length method. Transition rates matrix along a user defined direction and proportions of categories are computed.
tpfit_ml(data, coords, direction, tolerance = pi/8, mle = "avg")tpfit_ml(data, coords, direction, tolerance = pi/8, mle = "avg")
data |
a categorical data vector of length |
coords |
an |
direction |
a |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
mle |
a character value to pass to the function |
A 1-D continuous-lag spatial Markov chain is probabilistic model which involves a transition rate matrix computed for the direction . It defines the transition probability through the entry of the following matrix
where is a positive lag value.
To calculate entries of the transition rate matrix, we need to compute the mean lengths and the embedded transition probabilities.
By the use of the mean lengths, diagonal entries of are computed as
where is the mean length of the -th category.
The off-diagonal transition rates of the matrix are estimated by the use of embedded transition probabilities and mean lengths:
where is a specific embedded transition probability.
When some entries of the matrix are not identifiable, it is suggested to vary the tolerance coefficient or to set the input argument mle to "mlk".
An object of the class tpfit is returned. The function print.tpfit is used to print the fitted model. The object is a list with the following components:
coefficients |
the transition rates matrix computed for the user defined direction. |
prop |
a vector containing the proportions of each observed category. |
tolerance |
a numerical value which denotes the tolerance angle (in radians). |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.tpfit, print.tpfit, multi_tpfit_ml, transiogram
data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_ml(ACM$MAT5, ACM[, 1:3], c(0, 0, 1))data(ACM) # Estimate the parameters of a # one-dimensional MC model tpfit_ml(ACM$MAT5, ACM[, 1:3], c(0, 0, 1))
The function estimates transition probabilities matrices for a -D continuous lag spatial Markov chain.
transiogram(data, coords, direction, max.dist = Inf, mpoints = 20, tolerance = pi / 8, reverse = FALSE)transiogram(data, coords, direction, max.dist = Inf, mpoints = 20, tolerance = pi / 8, reverse = FALSE)
data |
a categorical data vector of length |
coords |
an |
direction |
a |
max.dist |
a numerical value which defines the maximum lag value. It's |
mpoints |
a numerical value which defines the number of lag intervals. |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
reverse |
a logical value. If |
Empirical probabilities are estimated by counting such pairs of observations which satisfy some properties, and by normalizing the result.
A generic pair of sample points and , where , must satisfy the following properties:
where is a non negative real value, while denotes the maximum lag value (max.dist) and is the number of lag intervals (mpoints).
the lag vector must have the same direction of the vector (direction) with a certain angular tolerance.
An object of the class transiogram is returned. The function print.transiogram is used to print computed probabilities. The object is a list with the following components:
Tmat |
a 3-D array containing the probabilities. |
LOSE |
a 3-D array containing the standard error calculated for the log odds of the transition probabilities. |
lags |
a vector containing one-dimensional lags. |
type |
a character string which specifies that computed probabilities are empirical. |
Luca Sartore [email protected]
Carle, S. F., Fogg, G. E. (1997) Modelling Spatial Variability with One and Multidimensional Continuous-Lag Markov Chains. Mathematical Geology, 29(7), 891-918.
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
predict.tpfit, predict.tpfit, plot.transiogram
data(ACM) # Estimate empirical transition # probabilities by points transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 200, 5)data(ACM) # Estimate empirical transition # probabilities by points transiogram(ACM$MAT3, ACM[, 1:3], c(0, 0, 1), 200, 5)
The function classifies points which appertain to a same directional line.
which_lines(coords, direction, tolerance = pi / 8)which_lines(coords, direction, tolerance = pi / 8)
coords |
an |
direction |
a |
tolerance |
a numerical value for the tolerance angle (in radians). It's |
The algorithm used by this function searches the nearest points to a directional line. The function classifies such pairs of points that have the minimum distance and the same direction of the vector .
This operation is done to order points, so that it's possible to compute mean lengths (mlen) and embedded transition probabilities (embed_MC).
A numerical vector containing the line number for each point.
Luca Sartore [email protected]
Sartore, L. (2010) Geostatistical models for 3-D data. M.Phil. thesis, Ca' Foscari University of Venice.
data(ACM) direction <- c(0,0,1) loc.id <- which_lines(ACM[, 1:3], direction)data(ACM) direction <- c(0,0,1) loc.id <- which_lines(ACM[, 1:3], direction)