解析の詳細

The Metropolis-Hastings Algorithm

In this program the particular algorithm used is the Metropolis-Hastings Algorithm.

The way this works is if you have a particular solution t^(j) and you wish to find the next solution t^(j+1) then you choose a trial move to t' and calculate the ratio:

a = |∂t'/∂t^(j)| p(t'|y) / p(t^(j)|y)

If a > 1 then t^(j+1) = t'
else with probability a set t^(j+1) = t'
and otherwise set t^(j+1) = t^(j)

In many cases of simple move the Jacobean |∂t'/∂t^(j)| is 1, but in some multiple-parameter moves it is not.

MCMC trial moves

There are three types of trial move used in the MCMC analysis. The first of these is where a single parameter is updated (order of parameter updating is randomised). Such a parameter will be given a trial value which lies somewhere in the range of possibilities provided by the constraints. To present these moves visually we will use a "." to represent a normal parameter and a "|" to denote a parameter which is a boundary of a group. The following series of representations shows single parameter updating in progress.

    |   .  .  ..  .   .   |  . .  .|  .    .    .  |
    |   ..    ..  .   .   |  . .  .|  .    .    .  |
    |   ..    ..  .   .   |  .   ..|  .    .    .  |
    |   ..    ..  .   .   |  .   ..  |.    .    .  |
    

The second form of move is where a group of events is moved together as in:

    |   ..    ..  .   .   |  .   ..  |.    .    .  |
    |   ..    ..  .   .|  .   ..  |   .    .    .  |
    |   ..    ..  .   .  |  .   ..  | .    .    .  |
    |   ..    ..  .   . |  .   ..  |  .    .    .  |
    

The final form of move is where a group of events is expanded or contracted by moving one of the boundaries as in:

    |   ..    ..  .   . |  .   ..  |  .    .    .  |
    |   ..    ..  .   . |  .   ..  |...|
    |   ..    ..  .   . |  .   ..  | .  .  . |
    |   ..    ..  .   . |  .   ..  |  .   .   .  |
    

These latter two moves are designed to improve mixing between states and convergence. These two moves are also sometimes combined where one boundary affects more than one group of events.

The last of these moves has a Jacobean which is not 1. If t_a and t_b are the parameters for the two boundaries in question and there are n events being expanded or contracted then:

|∂t'/∂t^(j)| = [(t_b'-t_a')/(t_b^(j)-_a^(j))]ⁿ

Passes and burn-in

The program is set up to perform at least 10 groups of MCMC trials. In each of these groups we:

• Initialise the variables (the order of initialisation is randomised)
• Sort the variables to find a solution compatible with the constraints
• Burn in the MCMC (10% of the total group length)
• Perform the main MCMC trials
• Save the distributions

The initial default settings are for a total number of 30,000 passes, that is 3,000 per group; each pass includes a trial of each variable and some multiple variable trials as discussed in the previous section. By comparing the results from each group of trials we decide if the trial number is high enough by looking at the convergence integrals (see below). If any of these are lower than 95% of the expected maximum, the number of passes is automatically increased by a factor of two.

The agreement indices are based on the ratio of the likelihood under two different models.

In the most simple model is where each parameter has a uniform prior, or the prior is independent of t. In such cases we have:

p(t|y) ∝ p(y|t)

Individual agreement indices

Often an observation y_i relates to a single parameter t_i and under this simple model as there are no relationships between the parameters we have:

p(t_i|y_i) ∝ p(y_i|t_i)

Under our full model t_i depends on all of the other parameters and observations through the equation:

p(t|y) ∝ p(y|t) p(t)

From this we can integrate over all other parameters to get the marginal probability density, p(t_i|y).

We then define the agreement index to be:

F_i = ∫p(y_i|t_i) p(t_i|y) dt_i / ∫p(y_i|t_i) p(y_i|t_i) dt_i
A_i = 100 F_i

If the posterior and likelihood distributions are the same this ratio is 1 (or 100%), in some instances it can also rise to greater than 1 (> 100%) where the posterior picks up the higher portions of the likelihood distribution.

In outlier analysis we have an associated parameter φ_i which is 1 when the measurement is an outlier and 0 when it is not. The prior for φ_i being 1 is q_i and for it to be zero is (1-q_i). When outlier analysis is in use the following factor is added into the value of F_i:

(q_i φ_i + (1-q_i)(1-φ_i))/(q_i² + (1-q_i)²)

This factor is on average one if the posterior for φ_i=1 is q_i but if the posterior outlier probablity is higher than the prior the factor is usually smaller (as long as q_i < 0.5). This modification of the agreement index (introduced in v4.1.4 results in lower agreement indices where there are more outliers which helps with model comparison. This factor makes no difference when outlier analysis is not used and in general if outliers are accepted as part of the model, the agreement indices are not needed (except for parameters which do not have outliers applied.

The individual agreement indices are reported in the A column of the results table.

Overall measures of agreement

We can expand on the treatment of single parameters, to look at the model as a whole. One way we can do this is to calculate the quantity:

F_model = ∫p(y|t) p(t|y) dt / ∫p(y|t) p(y|t) dt

Which is a kind of pseudo Bayes factor. Alternatively we can calculate the product of the individual agreement indices (which generally gives approximately the same value unless the parameters are very highly correlated):

F_overall = ∏_i A_i

These ratios tend to get further and further from 1 for larger number of observations. For this reason we define two indices based on these ratios which depend on the number of likelihood distributions n:

A_model = 100 F_model^1/√n
A_overall = 100 F_overall^1/√n

In principle F_model and A_model are better measures. A_overall is given for compatibility with previous versions of the program.

Combinations and choice of threshold

In the case of simple combinations generated from the function Combine(), D_Sequence or the & operator, there is only one independent parameter. In such cases F_overall will be equal to F_model. Because in this case there is only one independent parameter in the comparison model, we give the agreement index in this case a special name:

A_comb = 100 F_overall^1/√n

If we combine normal distributions in this way we can also do a χ² test and it turns out that the level at which this fails (95% null hypothesis) is when:

A_n = 100/√(2n)

We can also calculate the logarithmic average of the individual agreement indices that are required to fail the same test by calculating A'_n where:

A'_n = 100 (A_overall/100)^1/√n = 100 F_overall^1/n

If we do this for a range of values of n:

  ________________________

    n     An(%)    A'n(%)
  ________________________
    1     70.7     70.7
    2     50.0     61.3 
    3     40.8     59.6
    4     35.4     59.5
    5     31.6     59.8
    6     28.9     60.2
    7     26.7     60.7
    8     25.0     61.3
    9     23.6     61.8
   10     22.4     62.3
   15     18.3     64.5
   20     15.8     66.2
   25     14.1     67.6
   30     12.9     68.8
   40     11.2     70.7
   50     10.0     72.2
   60      9.1     73.4
   80      7.9     75.3
  100      7.1     76.7
 ________________________
    

the value of A'_n is around 60% for a wide range of values on n. This is the justification for taking 60% as the threshold for acceptability of individual agreement indices. The argument for the same threshold for A_overall and A_model is that we expect ln(F) to randomly walk from 0 as we increase the number of parameters. This threshold is given the symbol:

A'_c = 60%

The choice of threshold of acceptability for these indices is somewhat arbitrary - though no more so than the recommendations for Bayes factors. In practice, however, the levels suggested here have been found with experience to be about right in identifying models that are inconsistent with the observations, and individual samples that are clearly aberrant in their context.

The overlap integral used to test for convergence is calculated in this way. The marginal posterior density distribution is built up in a series of MCMC trial groups, each starting from a different initial state. If p_i(t_i|y) is the density accumulated in all groups of trials so far and p'_i(t_i|y) is the density function estimated from a new group of trials, the quantity:

C_i = 100 [∫ p'_i(t_i|y) p_i(t_i|y) dt_i]²/[(∫ p'_i²(t_i|y) dt_i)(∫ p_i²(t_i|y) dt_i)]

is used as an estimate of convergence. This should be equal to 100% if both distributions are identical.

Command line options

Options can also be set by command line options when calling the analysis program. Furthermore default options for an installation of the program are set in OxCal.dat file found in the same directory as the program files. The possible options are given below:

-afilename   append log to a file*
-b1          BP                     -b0     BC/AD
-cfilename   use calibration data file
-d1          plot distributions     -d0     no plot 
-fn          default iterations for sampling in thousands
-g1          +/-                    -g0     BC/AD/BP 
-h1          whole ranges           -h0     split ranges 
-in          resolution of n
-ln      limit on number of data points in calibration curve (see Resolution)
-m1          macro language         -m0     simplified entry
-n1          round ranges           -n0     no rounding
-o1          include converg info   -o0     do not include 
-p1          probability method     -p0     intercept method 
-q1          cubic interpolation    -q0     linear interpolation 
-rfilename   read input from a file+
-s11         1 sigma  ranges        -s10    range not found 
-s21         2 sigma ranges         -s20    range not found
-s31         3 sigma ranges         -s30    range not found
-t1          terse mode             -t0     full prompts
-u1          uniform span prior     -u0     as in OxCal v2.18 and previous
-v1          reverse sequence order -v0     chronological order
-wfilename   write log to a file*
-yn          round by n years       -y0     automatic rounding

    

* Note that with either of these options the tabbed results will then be sent to the console output and can therefore be redirected to a file or a pipe; the standard redirection > or > > can be used instead if only the log file needs redirecting.

+ Note that the standard redirection < can also be used.

Automatic range rounding

If automatic range rounding is selected, the rounding is carried out depending on the overall range of the result. The following table shows the degree of rounding selected.

______________________________________________

Total range        Round to the nearest
______________________________________________

   1 -   50          1 year
  50 -  100          5 years
 100 -  500         10 years
 500 - 1000         50 years
1000 - 5000        100 years
...                ...
______________________________________________        

 BP(AD(1800))+sigmaBP(AD(1800))*randN()

/* Example of a macro for calibrating all of
the radiocarbon dates between 1000 and 2000 in
increments of 10 */

// define and initialise the variable
var(a);a=1000; 
while(a<=2000) 
{ 
 // calibrate the date
 R_Date(a,30); 
 a=a+10; 
};