All distributions and calibration curves are stored at the resolution set in the system options, rs. There is also a calculation resolution defined which is rc=1 for rs=1...19, rc=10 for rs=20...199 and so on. All dates are rounded to this value (input and output). Interpolation between the stored points is linear.

When integrations or differentiations are carried out they are at the resolution rc. The details of the interpolation methods (such as methods of rounding used) have been carefully chosen to give the expected results and variation from the analytical values are rarely more than a single year with the standard options.

The files for the calibration curve usually have a different resolution to the internal storage resolution and so some form of interpolation is needed. This can either be linear or a cubic function depending on the setting in the system options. The cubic interpolation does not fit a spline function as this is very time consuming to calculate and can have some undesirable features such as large excursions between points. The cubic function used here gives a smooth curve with a continuous first differential but gives very little overall difference from the linear interpolation. The form of the function between two points is simply defined by the four surrounding points. If fj defines the function at tj the interpolation between tj and tj+1 is given by f(t) where:

The calibration curve is stored in two arrays one ri defining the radiocarbon age of the tree rings and another sigmai defining the errors associated with these measurements. Both ri and sigmai which are stored at the resolution rs are generated from the supplied calibration curves using the above interpolation method.

Calendar and BP dates

Output and input can be given in terms of calendar years yCAL (AD/BC) or years before 1950 yBP (cal BP). The relationship between these is simply:

yCAL = 1950 - yBP


10BP = 1940AD, 11950BP = 10000BC

It should be noted that this does imply a year 0 in the AD/BC sequence which is strictly speaking incorrect. With radiocarbon dates the problem is clearly semantic, with historical evidence it should be borne in mind that age differences across the BC/AD boundary are actually one year larger that they should be. Alternatively negative numbers (BC) should always be taken as the start of the year and positive numbers (AD) as the end. Thus -1 is the start of the first year BC whereas +1 is the end of the year 1AD. The reason for this problem is that in order to keep the internal representation of the numbers consistent it is very difficult to have to deal with a number set which goes from -1 to 1.

Radiocarbon calibration

The radiocarbon calibration itself is performed (using the verb R_Date) by a comparison of the measured radiocarbon age to the values stored in the curve (the methods used are similar to those used by Stuiver and Reimer 1993 and van der Plicht 1993; also see Dekling and van der Plicht 1993; the error terms in the calibration curve are taken into account as in the latest versions of CALIB). This allows a variance distribution to be calculated: if the radiocarbon date is rm with a measurement error sigmam the variance distribution vi and the probability distribution pi are then given by:

In this program the distribution is left normalised to a maximum of 1 rather than the actual probability of any individual year.

NOTE This is different to old versions of OxCal (pre 3.2) where pi was simply set to exp(-vi/2).

See also [Archaeological Considerations]

Reservoir corrections

Reservoir corrections can be made using the command Reservoir. This is suitable for a simple reservoir with a uniform diffusion rate with the atmosphere. Such a diffusion rate can be expressed in terms of a time constant tau. If the reservoir concentration is r(t) and the atmospheric concentration is given by R(t) we can write down the differential equation:

Solution of this differential equation requires a knowledge of the curve R(t) for all times before t. A linear extrapolation is assumed before the start of the curve using a gradient estimated from the first half of the curve R(t). The uncertainties in this are assumed to be ten times larger than those quoted for the first point in the curve; in practice these assumptions are unlikely to be significant unless the time constant is very long or you are considering points close to the start of the calibration curve.

Treatment of the uncertainties is more complicated. If the uncertainties associated with each point on the calibration curve are assumed to be independent the uncertainties in the reservoir curve should be smaller. In practice the errors almost certainly to some extent systematic. They have therefore been treated in exactly the same way as the concentrations themselves: if sigma(t) and Sigma(t) are the respective uncertainties we assume:

If there are also uncertainties in tau the solution of the equations would involve a double integration which would in practice be very slow. Another algorithm has therefore been adopted which is to increase the sigma in proportion to the difference between R(t) and r(t). Thus if the uncertainty in tau is deltatau:

For the oceans a properly modeled ocean curve should be used (see Stuiver et al 1998 - marine data). Local corrections can then be made using a Delta_R correction term:

See also [Calibration Data]

Mixed calibration curves

Calibration curves can be mixed using the Mix_Curves statement. This defines the two curves to be mixed and the proportion of the second one to be included. The mixing is calculated in terms of radiocarbon concentration rather than radiocarbon age. If the radiocarbon concentrations of the two curves are:

R1 E1
R2 E2

and the proportion of the second curve is P D then the resultant distribution is given by:

Rr = (1 - P) R1 + P R2
Er = Sqrt[((1 - P) E1)2 + (P E2)2 + (D (R2 - R1))2]

See also [Calibration Data]

Calendar dates and Asymmetric dates

Calendar dates are either entered as tcdtc with the command CAL or tc + d(+)tc - d(-)tc. In some cases such dates may be entered as accumulated doses dm from luminescence measurements in which case these are converted to calendar dates on input using the measurement year tm and dose rate dr and the simple relationships:

tc = tm - (dm/dr)

dtc = -(ddm/dr)

Rounding to the nearest rc will take place at this stage so you may notice a slight change in the entered values especially if rc is 10 or 100.

Calculation of symmetric probability distributions is simple:

The function used for asymmetric dates is rather more complex:

Proportional errors and factors

Proportional errors (Error) are particularly relevant to luminescence dating methods and proportional factors (Factor) can also be used for this type of application. Both of these are related to the year of measurement tm. A proportional factor f can be applied to a distribution with the mapping:

p'(t) = p(t/f)

And a proportional error df by using the mapping:

The distribution is then renormalised.

In the program these error factors are normally calculated before each distribution is reported except in the case of functions such as Combine which give a resultant distribution when the factor is only applied to the final result to prevent the systematic errors being reduced in the combination process.

Range calculation

Ranges are calculated to the resolution rc (that is normally to the nearest year) by linear interpolation of the probability or variance arrays. The boundaries of the age ranges for the intercept method are simply given by variance distribution: the levels required being 1, 4 and 9 for the 1, 2 and 3 standard deviation ranges (this is the method first used by Stuiver and Reimer 1986; see also Bowman 1990). This method is only used for radiocarbon dates and calendar dates, the probability method always being used for the results of any more complex analysis.

The probability method (selected for all types of distribution in the system options) calculates the ranges in a different way (similar to the method used by van der Plicht 1993). The elements of the probability distribution array pi are sorted by size and the integral normalised to 1.0. Starting from the top the array is then integrated until a certain proportion of the total is achieved (68.2%, 95.4% or 99.7%) and the level at this point in the distribution found pr. The ranges can then be defined as those parts of pi > pr.

If whole ranges are selected from the system options with the probability method a slightly different method is employed in order to generate floruits (see Aitchison et al 1991): the probability distribution is normalised to an integral of 1.0 and then the distribution is integrated from each end until a certain proportion of the curve has been excluded (15.9%, 2.3% or 0.15% from each end); the range defined is then the part of the distribution between these two points.

Integrated distributions (generated by the functions Before and After) define ranges directly from the height of the distribution using the values 0.682, 0.954 and 0.997.

Combinations and Wiggle Matching

Combinations of radiocarbon dates prior to calibration (R_Combine) and direct combination of calendar dates with gaussian errors (C_Combine) is performed in the normal way with a chi squared test (see for example Shennan 1988 p65).

Combinations of probability distributions (Combine) are simply done by using the Bayesian rules for combinations of probabilities (see Bayes 1763 and Doran and Hodgson 1975): if we have two probability distributions p1(t) and p2(t) these are combined as:

r(t) = p1(t) p2(t)

or more generally:

For the purposes of display the maximum of the resultant distributions is always normalised to 1.

If within a group defined for the function Combine the distributions are given a Gap gi then the combination is performed as:

We can then define a new set of original distributions p'i using

p'i = r(t+gi)

A very similar method to this is used for wiggle matching using the command D_Sequence, only difference being the way in which the gap is defined (between each successive distribution). A probability distribution r(t) is always calculated for the start of the sequence. This is given by:

and the resultant distributions then calculated using:

In the case of Bayesian wiggle matches and combinations, the program also calculates the chi-squared value for the best fit (ie the highest point on the probability distribution). This is reported in the text log file. For wiggle matching tree ring sequences, where the overall precision can be very high, you use a resolution of one year.

See also [Archaeological Considerations]

First and Last Dated Events in a Group

The probability of being after a single event is given by:

and so if a group of events are independent the probability of being after all of them is given by:

This is the distribution (normalised to a maximum of 1) returned by After. From this a distribution r'(t) can be calculated which gives a probability distribution for the last of the group of events:

r'(t) = dr(t)/dt

This is the distribution returned by the request Last within a phase if MCMC sampling is not needed.

The probabilities of being before a group of events and a distribution for the first event of a phase can be similarly defined.

WARNING: These methods assume that the events are entirely independant; in most cases a much better estimate will be arrived at using MCMC sampling from a phase which is enclosed within Boundary events.

See also [Archaeological Considerations]

Offset dates and Age Differences

If a distribution is offset by an amount dt the probability distribution is simply given by:

r(t) = p(t-dt)

If the offset has an error associated with it then dt+-sigma the distribution is given by:

A similar method is used to calculate a probability distribution for the age difference between two independent distributions (only employed when MCMC sampling is not necessary) using Difference.

And to shift one distribution by another using Shift:

See also [Archaeological Considerations]