OxCal > 解 析 > データの詳記

データの詳記



変数の型式

OxCalで使用できる変数型は,以下の3型式になります。

OxCalプログラムでは,いかなる解析過程においても定義した変数型を保持します。Number()とDate()は,変数型式を定義するために使用されるコマンドです。


時間軸の定義

複雑な解析や変数型に対応するため,基本時間軸の定義が非常に重要となります。自由度の高い時間軸という観点から,整数型の時間軸であるcal BPは,“年”のみの表現しかできないため基本軸としては不十分な暦です。OxCalプログラムの基本時間軸は,グレゴリオ暦を採用しています。

グレゴリオ暦では,BC/ADに直接関連させた連続する実数の時間軸(G)として表現することが可能です。ただし,BC/ADでは0年が存在しないため,G-1が2BCの1月1日00:00:00に対応しており,G1がAD1の1月1日00:00:00に対応しているので注意してください。

cal BPを時間軸として一般化するために,AD1950年の間にあたる日付を実際の0 cal BPとして定義しています(プログラムでは,G1950.5を使用)。

OxCalには,この様に暦を変換するための日付演算機能を備えてあります。(事前演算機能を参照)。以下に,その使用例を列挙してみました。

暦の詳細と暦の換算機能については,暦の定義を参照ください。

OxCalから出力されたプロットの横軸は,1年から開始するようになっています。表示するプロットの横軸にBC/ADの転換期が含まれる場合,“1BC/1AD”のように表現されます。世紀単位で出力されたプロットの横軸は以下のようになります。

     .     .     .     .     .     .     .
    300   200   100 1BC/1AD 101   201   301
    

年単位で表すと,横軸は以下のようになります。

     .     .     .     .     .     .     .
     3     2     1  1BC/1AD  2     3     4
    

上で挙げた例を実質的な見方をすると以下のような関係で表せます。

     .     .     .     .     .     .     .
     | 3BC | 2BC | 1BC | AD1 | AD2 | AD3 |     
    

上述した時間軸をグレゴリオ年で表現した場合は(±CEあるいはOxCalフォーマット),以下ように表せます。

     .     .     .     .     .     .     .
    -2    -1     0     1     2     3     4
    

また,BPで表現した場合は,各年の間が目盛に該当します。したがって,横軸は以下の通りです。

        .     .     .     .     .     .   
         1952  1951  1950  1949  1948  1947
    

BC/ADあるいはBCE/CE型式で年代あるいは年代期間を表示する場合,0年の存在しない暦でイベントの生じた年代を返します。ISO-8601(国際規格で定義され,天文学で採用)を選択している場合ではすべてのデータがCEで返され,0年も存在します。このとき,0年以前の年数は負でカウントされます。OxCalの基本軸であるグレゴリオ年を選択した場合では,浮動小数点を含んだ年代が返されます。すべての出力ファイルおよび各種データセットは,この型式が採用されています。

放射性炭素年代の較正のみにOxCalを使用している場合には,暦の定義について特別気を使う必要はありませんが,将来的には暦の定義で示した定義をもとに多様な年代型式で出力結果を返せるようにアップデートする予定です。


変数

各変数は,それぞれ変数名を割り当てて解析モデルに導入します。変数を定義するための主なコマンドは,以下の通りです。

Number()は,変数の型式で紹介したように数値型の実数を返すコマンドです。これは,次のセクションで扱うDate()と同様の機能を持ちます。Number以外のコマンドは,それぞれ確率密度関数を返します。関数を返すコマンドは,OxCalでは通常それ自体を尤度として利用されます。N()は正規分布,U()とTop_Hat()は一様分布,P()は任意の関数を扱え,Prior()はある事前情報から数値的に表現された分布を返します。

いずれの関数も,二通りの表現ができます。OxCalでは,Date()やNumber()のような変数型を宣言しなくても,自動的にモデルの内容から入力型式を認識します。下に挙げた例はすべて同じ意味になります。

    a=200;
    a=Number(200);
    Number("a",200);
    

等号を使った表現方法は,数学的に複雑な計算処理に向いています。ただし,複雑な変数名の設定はできませんし,+-/*()などの符号を挿入することもできません。変数型の定義例は以下の通りです。

// Simple numbers
n1=200;                  n2=166.66*1.2;

// Normal likelihood with a mean of 200 and a sigma of 20
a1=N(200,20);            N("a2",200,20); 

// Uniform likelihood between 180 and 200
b1=U(180,220);           U("b2",180,220); 
b3=Top_Hat(200,20);      Top_Hat("b4",200,20); 

// Exponentially falling likelihood with time constant 10
d1=P(0,100,exp(-d1/10)); P("d2",0,100,exp(-d2/10)); 

このように一般化した変数は,様々な解析で応用できます(確率密度分布の演算を参照)。

Maths ↓

The parameters of any Bayesian model are all treated in a similar way. Each parameter has a value ti and might also have observations or other information about it which are conceptually denoted as yi. Each parameter introduced is assumed to have a uniform uninformed prior p(ti). Where there is other information, this can be used to define a likelihood function p(yi|ti) for that parameter.

The functions N(), U(), Top_Hat(), P() and Prior(), can be used to directly assign such likelihoods:

Function Likelihood
N(ri,si) p(yi|ti) ∝ (1/(si√(2π))) exp(-(ti - ri)2/(2 si2))
U(ri,si) p(yi|ti) ∝ H(ti-ri)H(si-ti)
Top_Hat(ri,si) p(yi|ti) ∝ H(ti-ri+si)H(ri+si-ti)

where H(x) is the Heaviside step function which is 0 when x<0, 1/2 when x=0 and 1 when x>0. The P() function defines a likelihood function directly and the Prior() function can be used to provide such a function in numerical form. Usually the Prior() function is used to define a non-uniform prior but it can equally be used to provide a particular functional form for a likelihood where the prior is defined in some other way. Assuming no other prior information, a likelihood can be viewed as an informed prior, since for constant p(t):

p(ti|yi) ∝ p(yi|ti)p(ti) ∝ p(yi|ti)

日付

日付変数型に対するコマンドは,あるイベントの時間情報を定義し,日付の変数型を宣言します。主なコマンドは以下のとおりです。

Date()は,実数を日付として定義するためのコマンドです。実際の解析時には,OxCalフォーマットとして数値の変数型と同様に扱われます(暦の定義を参照)。変数で紹介したように尤度関数としても使用することができます。Age()もDate()と同様に実行できますが,ある年代よりも以前の日付を相対的に返すことができます。基準となる日付を設定していない場合,AD1950年の間を基準日付とします。Year=を使用すると,任意の数値を日付変数として関数や基準値に設定することができます。

R_Date(),R_Simulate(),R_Combine()は,次のセクションで扱います。

C_Date()は,OxCal ver.3から継承したもので,正規分布の尤度を返します。設定によりBPやグレゴリオ年(G)を使用して出力することもできます。C_Simulate()は,正規分布した尤度を得るための推定年代を返します。Sapwood()関数は,年輪年代のための辺材推定で扱います。

C_Combine()関数は,C_DateやC_Simulate()を組み合わせた年代を出力します。

下の例では,G1066.5年を様々なコマンドを利用して表現しています。

    a1=Date(1066.5);       a2=Date(AD(1066));        a3=Age(884);       a4=Age(934){Year=AD(2000);};
    b1=Date(N(1066.5,10)); b2=Date(N(AD(1066),10));  b3=Age(N(884,10)); b4=C_Date(AD(1066),10);
    c1=Date(U(AD(1066),AD(1093))); Date("c2",U(AD(1066),AD(1093)));
    d2=C_Simulate(AD(1066),10);    d3=Date(N(AD(1066)+randN()*10,10));

Maths ↓

Date parameters are treated in exactly the same way as any other parameters and mathematically the Date() and Number() functions are identical. The Age() function assumes that any value or likelihood distribution relates to the time before present (defined by the Year attribute), Y.

The following show the effect of the Age() function on the likelihood:

Function Likelihood
Age(N(ri,si)) p(yi|ti) ∝ (1/(si√(2π))) exp(-(Y-ti - ri)2/(2 si2))
Age(U(ri,si)) p(yi|ti) ∝ H(Y-ti-ri)H(si-Y+ti)
Age(Top_Hat(ri,si)) p(yi|ti) ∝ H(Y-ti-ri+si)H(ri+si-Y+ti)

The function C_Simulate() calculates a likelihood as:

Function Likelihood
C_Simulate(ri,si) p(yi|ti) ∝ (1/(si√(2π))) exp(-(ti - ri-ε)2/(2 si2))
where ε is sampled from N(0,si)

The function C_Combine() function evaluates a combined error weighted mean and standard error following the method of Ward and Wilson 1978. If we have several determinations ri with standard errors si we calculate a combined determination rc with standard error sc:

rc = (Σ ri/si2)/(Σ 1/si2)
sc = (Σ 1/si2)-1/2
T = Σ (ri-rc)2/si2

The likelihood distribution is then as for the C_Date() function. The T value has a χ2 distribution on n-1 degrees of freedom as discussed in Ward and Wilson 1978.

If there is an additional component of uncertainty sa, this is added as:

sc = √[(Σ 1/si2)-1+sa2]

放射性炭素年代の較正

較正年代の算出過程は,測定試料の物質や出土地点などによって様々なデータを較正や年代補正へ複合的に利用する必要があるため,他の年代測定法で行われているような年代値の補正と比較してはるかに複雑な構造をとっています。

OxCalは,このように複雑な放射性炭素のあらゆる年代較正に対応可能なプログラムを用意しています。較正に使用する主なコマンドは,以下の通りです。

較正年代

R_Date()は,OxCalにあらかじめ設定されている較正曲線に対し,通常の年代較正を行います。

R_F14C()は,Reimer et al. 2004で定義されているF14C値を測定値として入力するときに使用します。R_Simulateは,特定の年代を持つ試料に対して年代測定から得られるであろう測定値を返します。較正年代は,オプションからグレゴリオ年かcal BPの設定が変更できます。

R_Combine()は,複数データを組み合わせたり,新たに考慮が必要な系統誤差を付け加えることができます(例えば,単年性の植物を扱ったときには,8 14C年の誤差を加えるような場合。これについては,Stuiver et al. 1998)。また,R_Combine()では,カイ二乗検定も同時に行われます(Ward and Wilson 1978)。

  R_Date("OxA-2000",2000,20);
  R_Combine("Jar 1a")
  {
   R_Date("OxA-2001",2010,20);
   R_Date("OxA-2002",1970,20);
   R_Date("OxA-2003",2020,20);
  };
  R_Simulate("Mid AD10",AD(10),20);

これは,Bombシリーズの較正曲線を使用した例です。

 Options()
 {
  Resolution=0.2;
  Curve="Bomb04NH1.14c";
 };
 Plot()
 {
  Curve("Bomb04NH1","Bomb04NH1.14c")
  {
   Reservoir(0.25,0.1);
  };
  R_F14C("OxA-8000",1.345,0.003);
  R_Combine("Jar 1b")
  {
   R_F14C("OxA-8001",1.345,0.003);
   R_F14C("OxA-8002",1.347,0.003);
   R_F14C("OxA-8003",1.341,0.003);
  };
  R_Simulate("Mid 1980",AD(1980),20);
 };

Maths ↓

If the curve function is defined as a radiocarbon concentration r(t) ± s(t) then the calibrated likelihood distribution for a radiocarbon determination of ri ± si is proportional to:

p(yi|ti) ∝ exp[-(ri-r(ti))2/(2(si2+s2(ti)))]/√(si2+s2(ti))

In practice this calculation is performed at increments defined by the resolution (which is every five years by default). If we have a ΔR defined to be rd ± sd the the likelihood for the reservoir offset is given by:

p(yd|td) ∝ exp[-(rd-td)2/(2 sd2))]/sd

and the likelihood for the calibrated date ri ± si is given by:

p(yi|ti) ∝ exp[-(ri-r(ti)-td)2/(2(si2+s2(ti)))]/√(si2+s2(ti))

This is the likelihood that is used in OxCal during all MCMC analysis in line with the recommendations of Jones and Nicholls 2002. If the ΔR applies to only one calibration we can combine the two distributions to arrive at:

p(yi|ti,yd) ∝ exp[-(ri-r(ti)-rd)2/(2(si2+s2(ti)+sd2))]/√(si2+s2(ti)+sd2)

This is the distribution calculated in the calculation stage of the program, and if no MCMC analysis is required.

The function R_Simulate(ri,si) calculates a likelihood as:

p(yi|ti) ∝ exp[-(r(ri)-r(ti)-ε)2/(2(si2+s2(ti)))]/√(si2+s2(ti))
where ε is sampled from N(0,si)

The function R_Combine() function evaluates a combined error weighted mean and standard error following the method of Ward and Wilson 1978. If we have several determinations ri with standard errors si we calculate a combined determination rc with standard error sc:

rc = (Σ ri/si2)/(Σ 1/si2)
sc = (Σ 1/si2)-1/2
T = Σ (ri-rc)2/si2

The likelihood distribution is then as for the R_Date() function. The T value has a χ2 distribution on n-1 degrees of freedom as discussed in Ward and Wilson 1978.

If there is an additional component of uncertainty sa, this is added as:

sc = √[(Σ 1/si2)-1+sa2]

較正曲線

Curve()により使用する較正曲線を設定することができます。較正曲線の詳細については,較正曲線のセクションを参照していください。このコマンドは,較正曲線として,あるいは,複数データセットの比較のために使用することができます(比較するデータセットを用意している場合)。較正曲線の各データポイントは,三つの変数を持っています。

初期設定されている較正曲線の分解能はIntCal較正曲線(5年)の分解能と同じで,補間やデータの端数処理を行っていません。ただし5年よりも低い分解能へ設定し直した場合は,キュービック(または設定により直線)補間を行います。これは,データセットの比較においても同様です。設定した分解能より各データポイントの分解能が高い場合は,設定した解像度に合わせて処理します。

Delta_R()は,主に海洋リザーバを対象とし(海洋リザーバ較正曲線),年代較正前に測定値をΔRで補正します。ベイズ解析を目的とした解析モデルで使用すると,ΔR補正値が正規分布する変数として扱われます(関連する試料すべてに共通)。また,複数データの間にある誤差をDelta_R()で補正することもできます。例えばDelta_R()で宣言した変数をdとしたとき,放射性炭素量rmをR_Date(rm-d,sm)のように定義することで,dで補正した較正年代を求めることができます。

Mix_Curve()は,異なった炭素リザーバの混合割合とその誤差から較正曲線を組み合わせることができます。得られた較正曲線は,初期設定されているデータセットと同じように年代較正に利用できます。Mix_Curve()を含んだモデルでは,MCMCによる解析が常に行われています。ただし,第2較正曲線の寄与率は,0-100%以内で設定してください(混合割合は,通例10±10%のように設定します)。

基本的にモデル解析や年代較正は,放射性炭素濃度から計算が行われていることに注意してください。BPを優先して解析や年代較正を行いたい場合は,“F14Cの使用”をoffにしてください。この設定は,設定した解析や年代較正だけに適用され,混合リザーバなどの解析は放射性炭素濃度F14Cで実行されます。

年代較正に用いる較正曲線の選択は,較正曲線のダイアログ(単一較正のためには[Options > Curve],あるいは,入力ウィンドウの[Tools > Curves])から簡単に行えます。較正曲線に関するコマンドの使用例は以下の通りです。

    Curve("Atmospheric","IntCal09.14c");
    R_Date(2000,20);
    Curve("Oceanic","Marine09.14c");
    R_Date(2000,20);
    Delta_R("Local Marine",200,30);
    R_Date(2000,20);
    Curve("Decadal Turnover","IntCal09.14c")
    {
     Reservoir(10,5);
    };
    R_Date(2000,20);
    Mix_Curves("Mixed","Atmospheric","Local Marine",40,10);
    R_Date(2000,20);

Maths ↓

Curve construction

Where binning is required the method is as follows:

We need to take account of the fact theat there is often much real natural variability not reflected in the measurement uncertainty. We have the individual measurements rij with measurement uncertainties sij. The error weighted estimate of the mean is:

μ'i=(Σj rij/sij2)/(Σj (1/sij2))

And the average weighted variance (Bevington and Robinson 1992) is:

σi2= [(Σj rij2/sij2)/(Σj (1/sij2)) - μ'i2] [ni/(ni-1)]

The uncertainty of the mean is:

σμi2 = 1/(Σj (1/sij2))

So we take the variance si2 in the bin to be the larger of σμi2 and σi2. Overall for each bin we have three values which are:

ti=(Σj tij/sij2)/(Σj (1/sij2))
ri=(Σj rij/sij2)/(Σj (1/sij2))
si= √max([(Σj rij2/sij2)/(Σj (1/sij2)) - ri2]/[ni/(ni-1)] , 1/(Σj(1/sij2)))

The revised data-set is then interpolated onto the required points as described below.

If the data-points are more widely spread than the resolution requires, or they do not lie on the correct points, interpolation is performed. If linear interpolation is chosen (or if the interpolation is only for one point) then the function r(t) between two points (ti,ri,si) and (ti+1,ri+1,si+1) is given by:

ai = (ri+1 - ri)/(ti+1-ti)
r(t) = ti + ai(t-ti)

If cubic interpolation is used then the function is defined to pass through all data points and have a continuous value and first derivative. The interpolation is given by:

ai = (ri+1 - ri-1)/(ti+1-ti-1)
ai+1 = (ri+2 - ri)/(ti+2-ti)
δri = (ri+1 - ri)
δti = (ti+1 - ti)
ci = (3 δri - δti(2ai + ai+1))/δti2
di = ((ai+1 - ai) - (2 ci δti))/(3 δti2)
r(t) = ti + ai(t-ti) + bi(t-ti)2+ ci(t-ti)3

s(t) is interpolated in exactly the same way.

Curve mixing

The mixing ratio is treated as a parameter tmin Bayesian models. Where the mixing ratio is defined as rm ± sm this parameter is given a likelihood which is:

p(ym|tm) ∝ H(tm)H(1-tm)exp(-(tm-rm)2/(2sm2))

in other words it is a normal distribution truncated at 0 and 1. A posterior distribution for the mixing ratio will be generated which is based on the global model. The combined curve is given by:

r'(t) = (1-tm) r1(t) + tm r2(t)
s'(t) = (1-tm) s1(t) + tm s2(t)

Curve reservoir

If a reservoir time constant τ is set it is assumed that the reservoir contains carbon from a range of ages exponentially distributed with an average age of τ. If we assume that the atmospheric radiocarbon curve is given by r(t) and the reservoir curve by r'(t), we can see that:

r'(t) = (1/τ) ∫t r(u) exp(-(t-u)/τ) du

We also define the same relationship for the uncertainty:

s'(t) = (1/τ) ∫t s(u) exp(-(t-u)/τ) du

By numerical integration it is possible to calculate these functions. We assume a linear extrapolation of the curve to arrive at a starting point but with errors that are ten times higher than the first point in the curve. This approximation means that the resultant curve should not be relied upon within a few time constants of the start of the curve r(t).

Where there is an uncertainty in τ of δτ the additional uncertainty is calculated by numerically determining dr'(t)/dτ.

In order to give a representative result even when the errors are non-gaussian (as is the case with long time constants in the post-bomb period) we numerically evaluate the greater of:

dr'(t)/dτ ≈ [r'(t,τ+δτ)-r'(t,τ-δτ)]/(2δτ)
dr'(t)/dτ ≈ [r'(t,τ)-r'(t,τ-δτ)]/δτ
dr'(t)/dτ ≈ [r'(t,τ+δτ)-r'(t,τ)]/δτ

and r''(t) is found by averaging:

r''(t) ≈ [r'(t,τ+δτ) + r'(t,τ) + r'(t,τ-δτ)]/3

and adding the additional uncertainty in quadrature to ds'(t) to give a revised curve r''(t) and s''(t) such that:

r''(t)=r'(t)
s''(t)= √((s'(t))2 + (δτdr'(t)/dτ)2)

Note that this a slightly different way of calculating the uncertainty than in versions of OxCal prior to version 4. The differences are significant where the curve has rapid changes (such as with the post-1950 calibration curves).


年輪年代のための辺材推定

OxCalには,樹皮が欠落している樹木試料に対して,芯材/辺材境界から辺材年輪数(S)を推定するためのプログラムが実装されています。推定プログラムは,Dan Miles(Miles 2005)の研究論文をもとにしています。辺材推定に必要な変数は,以下の通りです。

ln(M)とln(R)に対して,ln(S)が直線近似関係にあります。つまり,芯材の平均年輪幅Mの増加と年輪数Rが正の相関を示し,辺材年輪数Sは平均年輪幅Mの増加によって負の相関を示します。辺材ツールでは,樹木の生育環境に合わせた解析プログラムにより直線近似を行います。

辺材推定の計算

OxCalでは,辺材年輪数が未知の樹木試料に対して辺材推定を行い年代を返すことができます(ただし,芯材/辺材の層が残存している場合に限ります)。辺材推定には,以下のコマンドが必要となります。

Sapwood_Model()では,辺材推定モデルに必要な変数を定義します。Sapwood()には,名前に続けて4つの変数をとります。

辺材推定プログラムでは,入力したデータをもとに尤度分布を返します。

Dan Milesは,スコットランドを含むイギリスを主体としたローマ帝国後をフィールドとした事例研究および環境データをもとに (Miles 2005),オークの年代測定を行っています。辺材推定コマンドの使用例は以下の通りです。

  Sapwood_Model("Mainland Britain", 2.77292, 0.100001, -0.275445,0.314286377);
  Sapwood("wa21", 1329, 243, 0, 1.06);
  Sapwood("wa22", 1354, 58, 6, 2.74);

辺材推定を行う場合,初期設定値の分解能(5年)では低すぎるため,[Tools > Options]から解像度を1年に設定してください。

Maths ↓

We define variables based on the logarithms of S, M and R:

s=ln(S), m=ln(M) and r=ln(R)

You can find a linear dependency of s on m and r. The residuals are normally distributed in s. This linear regression can be performed using the sapwood tool provided with this package, or any statistical package. The likelihood for S is given by:

p(yi|Si,a,br,bm) ∝ H(S) exp(-(a + brln(Ri) + bmln(Mi) - ln(Si))2/(2σ2))/Si

or if we define qi as the heartwood/sapwood boundary date for a sample, the likelihood for the felling date is given by;

p(yi|ti,a,br,bm) ∝ H(ti-qi) exp(-(a + brln(Ri) + bmln(Mi) - ln(ti-qi))2/(2σ2))/(ti-qi)

日付情報の変数化

放射性炭素年代や辺材推定では解析結果を尤度分布の関数型で出力するために,特別なコマンドをそれぞれ必要としてきました。しかし,この他の年代測定法では,最終的な年代測定値が正規分布をとることが少なくありません。このようなデータも,日付で紹介したコマンドで解析を行うことができます。

ただし,年代測定法によっては,得られた年代結果を正しくベイズ解析にかけられるように,不確定さを正確に表現しておく必要があります。ここでは事例として,ルミネッセンス年代測定法について紹介します。

今,4測定試料の推定蓄積線量と2つの年間線量率があるとします。ただし,測定試料における含水量の変動から蓄積線量が低く見積もられてしまうことがあります。

この場合,4つの測定試料の年代値をそれぞれ求めるために,7つの変数が必要となります。

OxCalへは以下のように変数を定義します。

  Year=2006;
  // define the independent parameters and set the likelihoods
  DE1=N(403.0,5);
  DE2=N(447.5,3);
  DE3=N(433.7,6);
  DE4=N(462.3,8);
  R1=N(0.08345,0.0009);
  R2=N(0.08467,0.0012);
  F=N(1.100,0.05);
  // calculate the dependent parameters
  D1=Age(DE1/(R1*F));
  D2=Age(DE2/(R1*F));
  D3=Age(DE3/(R2*F));
  D4=Age(DE4/(R2*F));

ここで示したように,一般的な変数から年代値を算出して編年研究に利用することができます。

全ての独立変数(DE1…DE4,R1,R2およびF)が,正規分布していることに注意してください。ただし,独立変数に従属する数値まで分布関数として扱うケースは,誤差が非常に大きくなる特別な場合を除いて,ほとんどありません。


相互参照

OxCalのver3までは,ある2つのグループに対して注目するイベントが許容か否かの適性を判断するためのXReferenceが提供されていました。ver4以降では,より一般的に相互参照を行うための方法が用意されています。

前のセクションまでは,解析モデル作成に向けた変数の宣言を中心に例を挙げてきました。解析を行っていく上では,特定のデータに対して相互参照が必要になる場合もあります。モデルの中にある2つのイベントが同時期であることを定義したり,二つの場所のΔR変数をそれぞれ繰り返し使用する場合は,以下のように解析します。

    Curve("Marine04","Marine04.14c");
    Delta_R("Region 1",200,30);
    R_Date("a",3000,30);
    Delta_R("Region 2",100,30);
    R_Date("b",3000,30);
    Delta_R("=Region 1");
    R_Date("c",3000,30);
    Delta_R("=Region 2");
    R_Date("d",3000,30);

相互参照では,参照したい変数の定義が重要になり,宣言した変数名から相互参照を行います。上の相互参照の例ではDelta_R()に同じ名前が複数使用されていますが,後半のコマンドには等号=を使って前に宣言した変数名を参照させています。このセクションの最後の例のように,&=演算子も使うこともできます。下の例では,開始年代が同時期かつ異なった終焉年代を持つ2つのグループの解析モデルです。

  Phase()
  {
   Sequence()
   {
    // first instance of the parameter "Start 1"
    Boundary("Start 1");
    Phase("1")
    {
     R_Date("a", 2000, 20);
     R_Date("b", 2010, 20);
     R_Date("c", 1920, 20);
     R_Date("d", 1900, 20);
     R_Date("e", 1903, 20);
    };
    Boundary("End 1");
   };
   Sequence()
   {
    // cross reference to the parameter "Start 1"
    Boundary("=Start 1");
    Phase("2")
    {
     R_Date("f", 2040, 20);
     R_Date("g", 2020, 20);
     R_Date("h", 1960, 20);
     R_Date("i", 2000, 20);
     R_Date("j", 1933, 20);
    };
    Boundary("End 2");
   };
  };

次の例は,上と同じデータをコマンドを使用せず表したものです(変数の宣言の仕方と相互参照のためのs1&=Boundary()に注目)。

  // define the likelihoods
  a=R_Date(2000, 20);
  b=R_Date(2010, 20);
  c=R_Date(1920, 20);
  d=R_Date(1900, 20);
  e=R_Date(1903, 20);
  f=R_Date(2040, 20);
  g=R_Date(2020, 20);
  h=R_Date(1960, 20);
  i=R_Date(2000, 20);
  j=R_Date(1933, 20);
  // define the model
  (s1=Boundary()) < (a | b | c | d | e) < (e1=Boundary());
  (s1&=Boundary()) < (f | g | h | i | j) < (e2=Boundary());

相互参照を行うことで,ひとつの独立変数から不確定さを正確に処理した解析モデルを作成することができます。また,宣言された変数に対して関連する従属変数を出力することもできます。