library(geostats)

Description

A list of documented functions may be viewed by typing help(package='geostats'). Detailed instructions are provided at https://github.com/pvermees/geostats/.

Author(s)

Maintainer: Pieter Vermeesch p.vermeesch@ucl.ac.uk

See Also

Useful links:

• https://github.com/pvermees/geostats/

A-CN-K compositions

Description

Synthetic A (Al_2O_3) – CN (CaO+Na_2O) – K (K_2O) data table.

Examples

data(ACNK,package='geostats')
ternary(ACNK,type='p',labels=c(expression('Al'[2]*'O'[3]),
                               expression('CaO+Na'[2]*'O'),
                               expression('K'[2]*'O')))

British coast

Description

A 512 \times 512 pixel image of the British coastline.

Examples

data(Britain,package='geostats')
p <- par(mfrow=c(1,2))
image(Britain)
fractaldim(Britain)
par(p)

rivers on Corsica

Description

A 512 \times 512 pixel image of the river network on Corsica.

Examples

data(Corsica,package='geostats')
p <- par(mfrow=c(1,2))
image(Corsica)
fractaldim(Corsica)
par(p)

detrital zircon U-Pb data

Description

Detrital zircon U-Pb data of 13 sand samples from China.

References

Vermeesch, P. “Multi-sample comparison of detrital age distributions.” Chemical Geology 341 (2013): 140-146.

Examples

data(DZ,package='geostats')
qqplot(DZ[['Y']],DZ[['5']])

A-F-M data

Description

FeO - (Na_2O + K_2O) - MgO compositions of 630 calc-alkali basalts from the Cascade Mountains and 474 tholeiitic basalts from Iceland. Arranged in F-A-M order instead of A-F-M for consistency with the ternary function.

Examples

data(FAM,package='geostats')
ternary(FAM[,-1])

Finnish lake data

Description

Table of 2327 Finnish lakes, extracted from a hydroLAKES database.

References

Lehner, B., and Doll, P. (2004), Development and validation of a global database of lakes, reservoirs and wetlands, Journal of Hydrology, 296(1), 1-22, doi: 10.1016/j.jhydrol.2004.03.028.

Examples

data(Finland,package='geostats')
sf <- sizefrequency(Finland$area)
size <- sf[,'size']
freq <- sf[,'frequency']
plot(size,freq,log='xy')
fit <- lm(log(freq) ~ log(size))
lines(size,exp(predict(fit)))

get the mode of a dataset

Description

Computes the most frequently occuring value in a sampling distribution.

Usage

Mode(x, categorical = FALSE)

Arguments

Value

a scalar

Examples

data(catchments,package='geostats')
m1 <- Mode(catchments$CaMg,categorical=TRUE)

m2 <- 1:50
for (i in m2){
   m2[i] <- Mode(rnorm(100),categorical=FALSE)
}
hist(m2)

Principal Component Analysis of 2D data

Description

Produces a 4-panel summary plot for two dimensional PCA for didactical purposes.

Usage

PCA2D(X)

Arguments

Examples

X <- rbind(c(-1,7),c(3,2),c(4,3))
colnames(X) <- c('a','b')
PCA2D(X)

calculate \bar{R}

Description

Given n circular or spherical measurements, the length of their normalised vector sum (\bar{R}) serves as a measure of directional concentration.

Usage

Rbar(trd, plg = 0, option = 0, degrees = FALSE)

Arguments

Value

a value between 0 and 1

Examples

data(striations,package='geostats')
Rbar(striations,degrees=TRUE)

\bar{R} to \kappa conversion

Description

Converts the empirical concentration parameter \bar{R} to the von-Mises concentration parameter \kappa.

Usage

Rbar2kappa(R, p = 1)

Arguments

Details

\bar{R} and \kappa are two types of concentration parameter that are commonly used in directional data analysis. \kappa is one of the parameters of the parametric von Mises distribution, which is difficult to estimate from the data. \bar{R} is easier to calculate from data. Rbar2kappa converts \bar{R} to \bar{\kappa} using the following approximate empirical formula:

\kappa = \frac{\bar{R}(p+1-\bar{R}^2)}{1-\bar{R}^2}

where p marks the number of parameters in the data space (1 for circle, 2 for a sphere).

Value

value(s) between 0 and +\infty

References

Banerjee, A., et al. “Clustering on the unit hypersphere using von Mises-Fisher distributions.” Journal of Machine Learning Research 6.Sep (2005): 1345-1382.

Examples

data(striations,package='geostats')
Rbar2kappa(Rbar(striations,degrees=TRUE))

additive logratio transformation

Description

Maps compositional data from an n-dimensional simplex to an (n-1)-dimensional Euclidean space with Aitchison's additive logratio transformation.

Usage

alr(dat, inverse = FALSE)

Arguments

Value

If inverse=FALSE, returns an (n-1) \times m matrix of logratios; otherwise returns an (n+1) \times m matrix of compositional data whose columns add up to 1.

Examples

xyz <- rbind(c(0.03,99.88,0.09),
             c(70.54,25.95,3.51),
             c(72.14,26.54,1.32))
colnames(xyz) <- c('a','b','c')
rownames(xyz) <- 1:3
uv <- alr(xyz)
XYZ <- alr(uv,inverse=TRUE)
xyz/XYZ

box counting

Description

Count the number of boxes needed to cover all the 1s in a matrix of 0s and 1s.

Count the number of boxes needed to cover all the 1s in a matrix of 0s and 1s.

Usage

boxcount(mat, size)

boxcount(mat, size)

Arguments

Value

an integer

an integer

Examples

g <- sierpinski(n=5)
boxcount(mat=g,size=16)
g <- sierpinski(n=5)
boxcount(mat=g,size=16)

Cantor set

Description

Calculates or plots a Cantor set of fractal lines, which is generated using a recursive algorithm that is built on a line segment whose middle third is removed. Each level of recursion replaces each black line by the same pattern.

Usage

cantor(n = 5, plot = FALSE, add = FALSE, Y = 0, lty = 1, col = "black", ...)

Arguments

Value

a square matrix with 0s and 1s.

Examples

plot(c(0,1),y=c(0,1),type='n',bty='n',ann=FALSE,xaxt='n',yaxt='n',xpd=NA)
cantor(n=0,Y=1.00,plot=TRUE,add=TRUE)
cantor(n=1,Y=0.75,plot=TRUE,add=TRUE)
cantor(n=2,Y=0.50,plot=TRUE,add=TRUE)
cantor(n=3,Y=0.25,plot=TRUE,add=TRUE)
cantor(n=4,Y=0.00,plot=TRUE,add=TRUE)

properties of 20 river catchments

Description

six different (three discrete, three continuous) measurements for twenty fictitious river catchments, containing their dominant lithology (categorical data), stratigraphic age (ordinal data), number of springs (count data), the pH of the river water (Cartesian quantity), its Ca/Mg ratio (Jeffreys quantity) and the percentage covered by vegetation (proportion).

Examples

data(catchments,package='geostats')
hist(catchments$pH)

plot circular data

Description

Plots directional data as ticks on a circle, with angles plotting in a clockwise direction from the top.

Usage

circle.plot(a, degrees = FALSE, tl = 0.1, ...)

Arguments

Value

no return value

Examples

data(striations,package='geostats')
circle.plot(striations,degrees=TRUE)

add points to a circular plot

Description

Adds directional data as points on an existing circle plot, with angles plotting in a clockwise direction from the top.

Usage

circle.points(a, degrees = FALSE, ...)

Arguments

Value

no return value

Examples

data(striations,package='geostats')
circle.plot(striations,degrees=TRUE)
md <- meanangle(striations,degrees=TRUE)
circle.points(md,pch=22,bg='black',degrees=TRUE)

centred logratio transformation

Description

Maps compositional data from an n-dimensional simplex to an n-dimensional Euclidean space with Aitchison's centred logratio transformation.

Usage

clr(dat, inverse = FALSE)

Arguments

Value

an n x m matrix

Examples

xyz <- rbind(c(0.03,99.88,0.09),
             c(70.54,25.95,3.51),
             c(72.14,26.54,1.32))
colnames(xyz) <- c('a','b','c')
rownames(xyz) <- 1:3
pc <- prcomp(clr(xyz))
biplot(pc)

colour plot

Description

Adds a colour bar to a scatter plot and/or filled contour plot. This function, which is based on base R's filled.contour function, is useful for visualising kriging results.

Usage

colourplot(
  x,
  y,
  z,
  X,
  Y,
  Z,
  nlevels = 20,
  colspec = hcl.colors,
  pch = 21,
  cex = 1,
  plot.title,
  plot.axes,
  key.title,
  key.axes,
  asp = NA,
  xaxs = "i",
  yaxs = "i",
  las = 1,
  axes = TRUE,
  frame.plot = axes,
  extra,
  ...
)

Arguments

Value

no return value

Examples

data('meuse',package='geostats')
colourplot(X=meuse$x,Y=meuse$y,Z=log(meuse$zinc),
           plot.title=title(main='Meuse',xlab='Easting',ylab='Northing'),
           key.title=title(main='log(Zn)'))

count the number of earthquakes per year

Description

Counts the number of earthquakes per year that fall within a certain time interval.

Usage

countQuakes(qdat, minmag, from, to)

Arguments

Value

a table with the number of earthquakes per year

Examples

data(declustered,package='geostats')
quakesperyear <- countQuakes(declustered,minmag=5.0,from=1917,to=2016)
table(quakesperyear)

declustered earthquake data

Description

Dataset of 28267 earthquakes between 1769 and 2016, with aftershocks and precursor events removed.

References

Mueller, C.S., 2019. Earthquake catalogs for the USGS national seismic hazard maps. Seismological Research Letters, 90(1), pp.251-261.

Examples

data(declustered,package='geostats')
quakesperyear <- countQuakes(declustered,minmag=5.0,from=1917,to=2016)
table(quakesperyear)

earthquake data

Description

Dataset of 20000 earthquakes between 2017 and 2000, downloaded from the USGS earthquake database (https://earthquake.usgs.gov/earthquakes/search/).

Examples

data(earthquakes,package='geostats')
gutenberg(earthquakes$mag)

ellipse

Description

Compute the x-y coordinates of an error ellipse.

Usage

ellipse(mean, cov, alpha = 0.05, n = 50)

Arguments

Value

a two-column matrix of plot coordinates

Examples

X <- rnorm(100,mean=100,sd=1)
Y <- rnorm(100,mean=100,sd=1)
Z <- rnorm(100,mean=100,sd=5)
dat <- cbind(X/Z,Y/Z)
plot(dat)
ell <- ellipse(mean=colMeans(dat),cov=cov(dat))
polygon(ell)

exponential transformation

Description

Map a logged kernel density estimate from [-\infty,+\infty] to [0,\infty] by taking exponents.

Usage

## S3 method for class 'density'
exp(x)

Arguments

Value

an object of class density

Examples

data(catchments,package='geostats')
lc <- log(catchments$CaMg)
ld <- density(lc)
d <- exp(ld)
plot(d)

fault orientation data

Description

Ten paired strike and dip measurements (in degrees), drawn from a von Mises - Fisher distribution with mean vector \mu=\{-1,-1,1\}/\sqrt{3} and concentration parameter \kappa=100.

Examples

data(fault,package='geostats')
stereonet(trd=fault$strike,plg=fault$dip,option=2,degrees=TRUE,show.grid=FALSE)

foram count data

Description

Planktic foraminifera counts in surface sediments in the Atlantic ocean.

Examples

data(forams,package='geostats')
abundant <- forams[,c('quinqueloba','pachyderma','incompta',
                      'glutinata','bulloides')]
other <- rowSums(forams[,c('uvula','scitula')])
dat <- cbind(abundant,other)
chisq.test(dat)

calculate the fractal dimension

Description

Performs box counting on a matrix of 0s and 1s.

Usage

fractaldim(mat, plot = TRUE, ...)

Arguments

Value

an object of class lm

Examples

g <- sierpinski(n=5)
fractaldim(g)

fractures

Description

A 512 \times 512 pixel image of a fracture network.

Examples

data(fractures,package='geostats')
p <- par(mfrow=c(1,2))
image(fractures)
fractaldim(fractures)
par(p)

create a Gutenberg-Richter plot

Description

Calculate a semi-log plot with earthquake magnitude on the horizontal axis,and the cumulative number of earthquakes exceeding any given magnitude on the vertical axis.

Usage

gutenberg(m, n = 10, ...)

Arguments

Value

the output of lm with earthquake magnitude as the independent variable (mag) and the logarithm (base 10) of the frequency as the dependent variable (lfreq).

Examples

data(declustered,package='geostats')
gutenberg(declustered$mag)

hills

Description

150 X-Y-Z values for a synthetic landscape that consists of three Gaussian mountains.

Examples

data(hills,package='geostats')
semivariogram(x=hills$X,y=hills$Y,z=hills$Z,model='gaussian')

Koch snowflake

Description

Calculates or plots a Koch set of fractal lines, which is generated using a recursive algorithm that is built on a triangular hat shaped line segment. Each level of recursion replaces each linear segment by the same pattern.

Usage

koch(n = 4, plot = TRUE, res = 512)

Arguments

Value

a res x res matrix with 0s and 1s

Examples

koch()

kriging

Description

Ordinary kriging interpolation of spatial data. Implements a simple version of ordinary kriging that uses all the data in a training set to predict the z-value of some test data, using a semivariogram model generated by the semivariogram function.

Usage

kriging(x, y, z, xi, yi, svm, grid = FALSE, err = FALSE)

Arguments

Value

either a vector (if grid=FALSE) or a matrix (if grid=TRUE) of kriging interpolations. In the latter case, values that are more than 10% out of the data range are given NA values.

Examples

data(meuse,package='geostats')
x <- meuse$x
y <- meuse$y
z <- log(meuse$cadmium)
svm <- semivariogram(x=x,y=y,z=z)
kriging(x=x,y=y,z=z,xi=179850,yi=331650,svm=svm,grid=TRUE)

Kolmogorov-Smirnov distance matrix

Description

Given a list of numerical vectors, fills a square matrix with Kolmogorov-Smirnov statistics.

Usage

ksdist(dat)

Arguments

Value

an object of class dist

Examples

data(DZ,package='geostats')
d <- ksdist(DZ)
mds <- cmdscale(d)
plot(mds,type='n')
text(mds,labels=names(DZ))

logistic transformation

Description

Maps numbers from [0,1] to [-\infty,+\infty] and back.

Usage

logit(x, ...)

## Default S3 method:
logit(x, inverse = FALSE, ...)

## S3 method for class 'density'
logit(x, inverse = TRUE, ...)

Arguments

Value

a vector with the same length of x

Examples

data(catchments,package='geostats')
lp <- logit(catchments$vegetation/100,inverse=FALSE)
ld <- density(lp)
d <- logit(ld,inverse=TRUE)
plot(d)

composition of Namib dune sand

Description

Major element compositions of 16 Namib sand samples.

References

Vermeesch, P. & Garzanti, E. “Making geological sense of ‘Big Data’ in sedimentary provenance analysis.” Chemical Geology 409 (2015): 20-27.

Examples

data(major,package='geostats')
comp <- clr(major)
pc <- prcomp(comp)
biplot(pc)

mean angle

Description

Computes the vector mean of a collection of circular measurements.

Usage

meanangle(trd, plg = 0, option = 0, degrees = FALSE, orientation = FALSE)

Arguments

Value

a scalar of 2-element vector with the mean orientation, either in radians (if degrees=FALSE), or in degrees.

Examples

data(striations,package='geostats')
meanangle(striations,degrees=TRUE)

Meuse river data set

Description

This data set gives locations and topsoil heavy metal concentrations, collected in a flood plain of the river Meuse, near the village of Stein (NL). Heavy metal concentrations are from composite samples of an area of approximately 15 m x 15 m. This version of the meuse dataset is a trimmed down version of the eponymous dataset from the sp dataset.

Examples

data(meuse,package='geostats')
semivariogram(x=meuse$x,y=meuse$y,z=log(meuse$cadmium))

palaeomagnetic data

Description

Ten paired magnetic declination (azimuth) and inclination (dip) measurements, drawn from a von Mises - Fisher distribution with mean vector \mu=\{2,2,1\}/3 and concentration parameter \kappa=200.

Examples

data(palaeomag,package='geostats')
stereonet(trd=palaeomag$decl,plg=palaeomag$incl,degrees=TRUE,show.grid=FALSE)

pebble orientations

Description

Orientations (in degrees) of 20 pebbles.

Examples

data(pebbles,package='geostats')
circle.plot(pebbles,degrees=TRUE)
m <- meanangle(pebbles,option=0,orientation=TRUE)
circle.points(m,degrees=TRUE,pch=22,bg='white')

3-magnet pendulum experiment

Description

Simulates the 3-magnet pendulum experiment, starting at a specified position with a given start velocity.

Usage

pendulum(
  startpos = c(-2, 2),
  startvel = c(0, 0),
  src = rbind(c(0, 0), c(0.5, sqrt(0.75)), c(1, 0)),
  plot = TRUE
)

Arguments

Value

the end position of the pendulum

Examples

p <- par(mfrow=c(1,2))
pendulum(startpos=c(2.1,2))
pendulum(startpos=c(1.9,2))
par(p)

generate bivariate random data

Description

Returns bivariate datasets from four synthetic distributions that have the shape of a circle, arrow, square and ellipse.

Usage

randy(pop = 1, n = 250)

Arguments

Value

a [2xn] matrix of random numbers

Examples

p <- par(mfrow=c(1,4))
for (i in 1:4){
   plot(randy(pop=i))
}
par(p)

Rb-Sr data

Description

Synthetic dataset of 8 Rb-Sr analysis that form a 1Ga isochron.

Examples

data(rbsr,package='geostats')
plot(rbsr[,'RbSr'],rbsr[,'SrSr'])
fit <- lm(SrSr ~ RbSr,data=rbsr)
abline(fit)

Spurious correlation

Description

Calculate the ‘null correlation’ of ratios, using the the spurious correlation formula of Pearson (1897).

Usage

rwyxz(
  mw,
  mx,
  my,
  mz,
  sw,
  sx,
  sy,
  sz,
  rwx = 0,
  rwy = 0,
  rwz = 0,
  rxy = 0,
  rxz = 0,
  ryz = 0
)

ryxy(mx, my, sx, sy, rxy = 0)

rxzyz(mx, my, mz, sx, sy, sz, rxy = 0, rxz = 0, ryz = 0)

Arguments

Value

the null correlation coefficient

References

Pearson, K. “Mathematical contributions to the theory of evolution. – on a form of spurious correlation which may arise when indices are used in the measurement of organs.” Proceedings of the Royal Society of London 60.359-367 (1897): 489-498.

Examples

rxzyz(mx=100,my=100,mz=100,sx=1,sy=1,sz=10)

semivariogram

Description

Plots the semivariance of spatial data against inter-sample distance, and fits a spherical equation to it.

Usage

semivariogram(
  x,
  y,
  z,
  bw = NULL,
  nb = 13,
  plot = TRUE,
  fit = TRUE,
  model = c("spherical", "exponential", "gaussian"),
  ...
)

Arguments

Value

returns a list with the estimated semivariances at different distances for the data, and (if fit=TRUE), a vector with the sill, nugget and range.

Examples

data(meuse,package='geostats')
semivariogram(x=meuse$x,y=meuse$y,z=log(meuse$cadmium))

Sierpinski carpet

Description

Returns a matrix of 0s and 1s that form a Sierpinski carpet. This is a two dimensional fractal, which is generated using a recursive algorithm that is built on a grid of eight black squares surrounding a white square. Each level of recursion replaces each black square by the same pattern.

Usage

sierpinski(n = 5)

Arguments

Value

a square matrix with 0s and 1s.

Examples

g <- sierpinski(n=5)
image(g,col=c('white','black'),axes=FALSE,asp=1)

calculate the size-frequency distribution of things

Description

Count the number of items exceeding a certain size.

Usage

sizefrequency(dat, n = 10, log = TRUE)

Arguments

Value

a data frame with two columns size and frequency

Examples

data(Finland,package='geostats')
sf <- sizefrequency(Finland$area)
plot(frequency~size,data=sf,log='xy')
fit <- lm(log(frequency) ~ log(size),data=sf)
lines(x=sf$size,y=exp(predict(fit)))

calculate the skewness of a dataset

Description

Compute the third moment of a sampling distribution.

Usage

skew(x)

Arguments

Value

a scalar

Examples

data(catchments,package='geostats')
skew(catchments$vegetation)

stereonet

Description

Plots directional data on a Wulff or Schmidt stereonet. The Wulff equal angle polar Lambert projection preserves the shape of objects and is often used to visualise structural data. The Schmidt equal area polar Lambert projection preserves the size of objects and is more popular in mineralogy.

Usage

stereonet(
  trd,
  plg,
  coneAngle = rep(10, length(trd)),
  option = 1,
  wulff = TRUE,
  add = FALSE,
  degrees = FALSE,
  show.grid = TRUE,
  grid.col = "grey50",
  tl = 0.05,
  type = "p",
  labels = 1:length(trd),
  pch = 21,
  bg = c("black", "white"),
  lty = c(1, 2),
  ...
)

Arguments

Author(s)

based on a MATLAB script written by Nestor Cardozo.

References

Allmendinger, R.W., Cardozo, N., and Fisher, D.M. “Structural geology algorithms: Vectors and tensors”. Cambridge University Press, 2011.

Examples

stereonet(trd=c(120,80),plg=c(10,30),degrees=TRUE,pch=16)
stereonet(trd=c(120,80),plg=c(10,30),degrees=TRUE,
          option=4,coneAngle=c(5,10),add=TRUE)

directions of glacial striations

Description

Directions (in degrees) of 30 glacial striation measurements from Madagascar.

Examples

data(striations,package='geostats')
circle.plot(striations,degrees=TRUE)

ternary diagrams

Description

Plot points, lines or text on a ternary diagram.

Usage

ternary(xyz = NULL, f = rep(1, 3), labels, add = FALSE, type = "p", ...)

Arguments

Examples

data(ACNK,package='geostats')
ternary(ACNK,type='p',labels=c(expression('Al'[2]*'O'[3]),
                               expression('CaO+Na'[2]*'O'),
                               expression('K'[2]*'O')))

composition of a further 147 oceanic basalts

Description

Major element compositions of 64 island arc basalts (IAB), 23 mid oceanic ridge basalts (MORB) and 60 ocean island basalts (OIB). This dataset can be used to test supervised learning algorithms.

References

Vermeesch, P. “Tectonic discrimination diagrams revisited.” Geochemistry, Geophysics, Geosystems 7.6 (2006).

Examples

library(MASS)
data(training,package='geostats')
ld <- lda(x=alr(training[,-1]),grouping=training[,1])
data(test,package='geostats')
pr <- predict(ld,newdata=alr(test[,-1]))
table(test$affinity,pr$class)

composition of 646 oceanic basalts

Description

Major element compositions of 227 island arc basalts (IAB), 221 mid oceanic ridge basalts (MORB) and 198 ocean island basalts (OIB). This dataset can be used to train supervised learning algorithms.

References

Vermeesch, P. “Tectonic discrimination diagrams revisited.” Geochemistry, Geophysics, Geosystems 7.6 (2006).

Examples

library(MASS)
data(training,package='geostats')
ld <- lda(x=alr(training[,-1]),grouping=training[,1])
pr <- predict(ld)
table(training$affinity,pr$class)

von Mises distribution

Description

Returns the probability density of a von Mises distribution, which describes probability distributions on a circle using the following density function:

\frac{\exp(\kappa\cos(x-\mu))}{2\pi I_0(\kappa)}

where I_0(\kappa) is a zero order Bessel function.

Usage

vonMises(a, mu = 0, kappa = 1, degrees = FALSE)

Arguments

Value

a scalar or vector of the same length as angles

Examples

plot(x=c(-1,1.2),y=c(-1,1.2),type='n',
     axes=FALSE,ann=FALSE,bty='n',asp=1)
a <- seq(from=-pi,to=pi,length.out=200)
d <- vonMises(a=a,mu=pi/4,kappa=5)
symbols(x=0,y=0,circles=1,add=TRUE,inches=FALSE,xpd=NA,fg='grey50')
lines(x=(1+d)*cos(a),y=(1+d)*sin(a),xpd=NA)

world population

Description

The world population from 1750 until 2014.

Examples

data(worldpop,package='geostats')
plot(worldpop)

get x,y plot coordinates of ternary data

Description

Helper function to generate bivariate plot coordinates for ternary data.

Usage

xyz2xy(xyz)

Arguments

Value

an n x 2 numerical matrix

Examples

xyz <- rbind(c(1,0,0),c(0,1,0),c(0,0,1),c(1,0,0))
xy <- xyz2xy(xyz)
plot(xy,type='l',bty='n')

Linear regression of X,Y-variables with correlated errors

Description

Implements the unified regression algorithm of York et al. (2004) which, although based on least squares, yields results that are consistent with maximum likelihood estimates of Titterington and Halliday (1979).

Usage

york(dat, alpha = 0.05, plot = TRUE, fill = NA, ...)

Arguments

Details

Given n pairs of (approximately) collinear measurements X_i and Y_i (for 1 \leq i \leq n), their uncertainties s[X_i] and s[Y_i], and their covariances cov[X_i,Y_i], the york function finds the best fitting straight line using the least-squares algorithm of York et al. (2004). This algorithm is modified from an earlier method developed by York (1968) to be consistent with the maximum likelihood approach of Titterington and Halliday (1979).

Value

A two-element list of vectors containing:

coef

the intercept and slope of the straight line fit

cov

the covariance matrix of the coefficients

References

Titterington, D.M. and Halliday, A.N., 1979. On the fitting of parallel isochrons and the method of maximum likelihood. Chemical Geology, 26(3), pp.183-195.

York, Derek, et al., 2004. Unified equations for the slope, intercept, and standard errors of the best straight line. American Journal of Physics 72.3, pp.367-375.

Examples

data(rbsr,package='geostats')
fit <- york(rbsr)