Extracting probability density function from specific depth changes model?
Zoë Thomas
Dear OxCal community
I’ve been struggling with this problem for several weeks now and I’m hoping that someone might be able to assist.
I have two sediment sequences, each with their own age model. The sediment sequences have a common volcanic ash layer. I’d like to extract the probability density function of a volcanic ash layer depth from each age model.
To do this, I simply add these 4 lines at the appropriate place in the code:
Date("Volcanic ash")
{
z=130;
};
In part, this works as intended, and allows me to extract the probability density function at depth 130.
However, when I compare the age model before and after I added this short piece of code, the age model looks different, and I believe that they are more different than would be expected due to the Bayesian process.
As an example, in one of the sediment sequences, I’m interested in looking at the probability density function from 5 volcanic ash layer depths. When the code to extract the pdf from 5 depths is added, this changes the agreement index of the model substantially (from over 100 to around 10). It seems that the short piece of code above is somehow changing the model. Should it do this? Is there an alternative way of extracting the probability density function at a specific depth?
I’m happy to share the code, but the age models take a while to run. I created a test age model (shorter, less dates) to try to recreate the problem, but the difference was far less pronounced.
Any help gratefully received.
All the best
Zoë Thomas
Christopher Ramsey
That does sound unusual and without the actual code it is hard to see why.
However, although the extra date command does not seem to change the model, it does in a subtle way which might be significant if the sedimentation rate is very variable.
This is ultimately due to the variation in the number of events/unit depth. To take an extreme example, if you put in one hundred events down a 100cm core the minimum k value would be 1 and the depth model would be fairly rigid. With only 4 depths it would allow it to be very flexible. This is ultimately because if you know there are this many events, in sequence, within the deposit we expect it to be more uniform in deposition rate.
So in your case, if the number of dates is very low, adding your five date significantly increases the total number of "events" and this might explain it. In this case the model with the layers included is really the correct one - rather than the one with only the dated samples, because you do know something happened at these points in the sequence.
I am more surprised it is affecting the agreement index to this extent - but presumably this is because the more rigid model imposed is in some way incompatible with the dates.
Best wishes
Christopher
> You received this message because you are subscribed to the Google Groups "OxCal" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to oxcal+un...@googlegroups.com.
> To view this discussion on the web visit https://groups.google.com/d/msgid/oxcal/5b5cc144-8b33-4e0d-89d3-8330a32e65bcn%40googlegroups.com.
Zoë Thomas
Dear Christopher
Thank you for your thoughtful response. I do see what you mean. However, I'm still quite surprised because the age model has quite a lot of dates (relatively speaking: 31 dates for 500cm). Note that these are not visible volcanic ash layers but cryptotephra, so I’m not expecting the sedimentation rate to change significantly, if at all. Essentially, I’d rather the age model didn’t change because I’d like to know how that particular depth is constrained by the radiocarbon age model surrounding it.
I have attached screenshots showing the plot vs. depth and have appended the code (with and without the ‘depth’ command). I have selected 3 volcanic ash layers to demonstrate (called H2, MB1, R1).
The agreement index (overall) for WITHOUT volcanic ash depth is 112.
The agreement index (overall) for WITH volcanic ash depth is 6.8.
I wondered about forcing it to think that it’s an outlier (and make up an age) but I don’t think that will give me the posterior probability distribution.
I’m stumped!
Many thanks for your help
Zoë
Code 1 (Without volcanic ash depths):
Options()
{
Curve="ShCal20.14c";
};
Plot()
{
Outlier_Model("General",T(5),U(0,4),"t");
P_Sequence("HP20",2,2,U(-2,2))
{
Boundary();
R_Date("Beta-193403",13630,140)
{
Outlier(0.05);
z=468;
};
Combine("HP20 29cm and Beta-241338")
{
Outlier("General", 0.05);
z=454.5;
R_Date("HP20 29cm", 13330,63);
R_Date("Beta-241338", 13320,50);
};
R_Date("HP20 19cm",12875,90)
{
Outlier("General", 0.05);
z=444.5;
};
R_Date("HP20 15cm",12393,72)
{
Outlier("General", 0.05);
z=440.5;
};
R_Date("HP20 12cm",12372,84)
{
Outlier("General", 0.05);
z=437.5;
};
R_Date("UBA-42351",11838,47)
{
color="red";
Outlier();
z=436;
};
R_Date("HP20 8cm",12213,69)
{
Outlier("General", 0.05);
z=433.5;
};
R_Date("HP20 4cm",12160,72)
{
Outlier("General", 0.05);
z=429.5;
};
R_Date("UBA-43463",11738,75)
{
Outlier(0.05);
z=427.5;
};
Combine("HP20 1cm")
{
Outlier("General", 0.05);
z=426.5;
R_Date("HP20 1cm", 11779,81);
R_Date("HP20 1cm", 11755,81);
};
R_Date("Beta-193402",10370,60)
{
Outlier(0.05);
z=400;
};
R_Date("Beta-241336",9940,40)
{
Outlier(0.05);
z=394.5;
};
R_Date("Beta-241335",10030,40)
{
Outlier(0.05);
z=390.5;
};
R_Date("Beta-193401",9250,80)
{
Outlier(0.05);
z=376;
};
R_Date("UBA-42350",8498,36)
{
Outlier(0.05);
z=369;
};
R_Date("UBA-43462",7587,38)
{
Outlier(0.05);
z=355;
};
R_Date("Beta-241334",7390,40)
{
Outlier(0.05);
z=342.5;
};
R_Date("UBA-42349",5796,29)
{
color="red";
Outlier();
z=338;
};
R_Date("UBA-42348",5953,29)
{
color="red";
Outlier();
z=335;
};
R_Date("UBA-43461",6029,37)
{
Outlier(0.05);
z=318;
};
R_Date("Beta-193400",5700,40)
{
Outlier(0.05);
z=301;
};
R_Date("UNSW-1197",4825,20)
{
Outlier(0.05);
z=285;
};
R_Date("UNSW-1196",4508,20)
{
Outlier(0.05);
z=260;
};
R_Date("UNSW-1195",4192,20)
{
Outlier(0.05);
z=243;
};
R_Date("UNSW-1194",3770,20)
{
Outlier(0.05);
z=220;
};
R_Date("UNSW-1193",3489,20)
{
Outlier(0.05);
z=197;
};
R_Date("UNSW-1192",2751,20)
{
Outlier(0.05);
z=149;
};
R_Date("UNSW-1191",1928,20)
{
Outlier(0.05);
z=95;
};
R_Date("UNSW-1190",1254,20)
{
Outlier(0.05);
z=65;
};
R_Date("UNSW-1189",1019,20)
{
Outlier(0.05);
z=55;
};
R_Date("UNSW-1188",444,20)
{
Outlier(0.05);
z=32;
};
Date("Top",N(2005))
{
z=0;
};
Boundary();
};
};
Code 2: (with volcanic ash depths):
Options()
{
Curve="ShCal20.14c";
};
Plot()
{
Outlier_Model("General",T(5),U(0,4),"t");
P_Sequence("HP20",2,2,U(-2,2))
{
Boundary();
R_Date("Beta-193403",13630,140)
{
Outlier(0.05);
z=468;
};
Combine("HP20 29cm and Beta-241338")
{
Outlier("General", 0.05);
z=454.5;
R_Date("HP20 29cm", 13330,63);
R_Date("Beta-241338", 13320,50);
};
R_Date("HP20 19cm",12875,90)
{
Outlier("General", 0.05);
z=444.5;
};
R_Date("HP20 15cm",12393,72)
{
Outlier("General", 0.05);
z=440.5;
};
R_Date("HP20 12cm",12372,84)
{
Outlier("General", 0.05);
z=437.5;
};
Date("HP20 - R1")
{
z=437.5;
};
R_Date("UBA-42351",11838,47)
{
color="red";
Outlier();
z=436;
};
R_Date("HP20 8cm",12213,69)
{
Outlier("General", 0.05);
z=433.5;
};
R_Date("HP20 4cm",12160,72)
{
Outlier("General", 0.05);
z=429.5;
};
R_Date("UBA-43463",11738,75)
{
Outlier(0.05);
z=427.5;
};
Combine("HP20 1cm")
{
Outlier("General", 0.05);
z=426.5;
R_Date("HP20 1cm", 11779,81);
R_Date("HP20 1cm", 11755,81);
};
R_Date("Beta-193402",10370,60)
{
Outlier(0.05);
z=400;
};
R_Date("Beta-241336",9940,40)
{
Outlier(0.05);
z=394.5;
};
R_Date("Beta-241335",10030,40)
{
Outlier(0.05);
z=390.5;
};
R_Date("Beta-193401",9250,80)
{
Outlier(0.05);
z=376;
};
R_Date("UBA-42350",8498,36)
{
Outlier(0.05);
z=369;
};
Date("HP05 - MB1")
{
z=367.5;
};
R_Date("UBA-43462",7587,38)
{
Outlier(0.05);
z=355;
};
R_Date("Beta-241334",7390,40)
{
Outlier(0.05);
z=342.5;
};
R_Date("UBA-42349",5796,29)
{
color="red";
Outlier();
z=338;
};
R_Date("UBA-42348",5953,29)
{
color="red";
Outlier();
z=335;
};
R_Date("UBA-43461",6029,37)
{
Outlier(0.05);
z=318;
};
R_Date("Beta-193400",5700,40)
{
Outlier(0.05);
z=301;
};
R_Date("UNSW-1197",4825,20)
{
Outlier(0.05);
z=285;
};
R_Date("UNSW-1196",4508,20)
{
Outlier(0.05);
z=260;
};
R_Date("UNSW-1195",4192,20)
{
Outlier(0.05);
z=243;
};
R_Date("UNSW-1194",3770,20)
{
Outlier(0.05);
z=220;
};
Date("HP05 - H2")
{
z=213.5;
};
R_Date("UNSW-1193",3489,20)
{
Outlier(0.05);
z=197;
};
R_Date("UNSW-1192",2751,20)
{
Outlier(0.05);
z=149;
};
R_Date("UNSW-1191",1928,20)
{
Outlier(0.05);
z=95;
};
R_Date("UNSW-1190",1254,20)
{
Outlier(0.05);
z=65;
};
R_Date("UNSW-1189",1019,20)
{
Outlier(0.05);
z=55;
};
R_Date("UNSW-1188",444,20)
{
Outlier(0.05);
z=32;
};
Date("Top",N(2005))
{
z=0;
};
Boundary();
};
};
Christopher Ramsey
Ok - it is indeed changing the k-profile. If you look at the distribution for k it is essentially uniform in the second case. The reason for this is that you have two elements at the same depth - so I think:
Options()
{
Curve="ShCal20.14c";
};
Plot()
{
Outlier_Model("General",T(5),U(0,4),"t");
P_Sequence("HP20",2,2,U(-2,2))
{
Boundary();
R_Date("Beta-193403",13630,140)
{
Outlier(0.05);
z=468;
};
Combine("HP20 29cm and Beta-241338")
{
Outlier("General", 0.05);
z=454.5;
R_Date("HP20 29cm", 13330,63);
R_Date("Beta-241338", 13320,50);
};
R_Date("HP20 19cm",12875,90)
{
Outlier("General", 0.05);
z=444.5;
};
R_Date("HP20 15cm",12393,72)
{
Outlier("General", 0.05);
z=440.5;
};
R_Date("HP20 12cm",12372,84)
{
Outlier("General", 0.05);
z=437.5;
};
{
z=437.5;
};*/
{
color="red";
Outlier();
z=436;
};
R_Date("HP20 8cm",12213,69)
{
Outlier("General", 0.05);
z=433.5;
};
R_Date("HP20 4cm",12160,72)
{
Outlier("General", 0.05);
z=429.5;
};
R_Date("UBA-43463",11738,75)
{
Outlier(0.05);
z=427.5;
};
Combine("HP20 1cm")
{
Outlier("General", 0.05);
z=426.5;
R_Date("HP20 1cm 14C2", 11755,81);
I've also renamed your combined radiocarbon dates - so that you don't get error messages about duplicate labels. You might wish to use R_Combine here in place of Combine.
Best wishes
Christopher
> <Zoe Code 1 OxCal without depth extract.png><Zoe Code 2 OxCal with depth extract.png>