Interpreting k from P_sequence

[Matthew Bolton]

Greetings,

I've been getting my feet wet with OxCal and have enjoyed the easy-to-use calibration and modelling.

However, I've got a bit of a quandary regarding the k value for a simple P-sequence deposition model. The profile I'm working on is an alluvial (flood deposit) sequence. Our knowledge of flood return intervals indicates that on average the k shouldn't be higher than 4.89 events per cm, and is probably considerably lower.

I have run models with variable k values ranging from my upper limit (4.89) to 0.5. However, my resultant k probability plots still indicate log(k) values less than zero. If I'm reading this right (and tell me if I'm not) my optimum k is actually a negative value?!? Attached is what appears to be the "best" result plot, from k0=0.5.

Further, I understand that so long as the model-selected k value is within the range allowed by the model, in my case, the standard, U(-2,2), the program should select the "right" k. However, I'd like to get an idea of the predicted event-rate to compare to stratigraphic analysis.

In short: a) How can I tell how low I have to make my k to yield a reliable model and b) if the model selects my k from a range, how can I tell what the resultant model "optimum" is?

The code I'm working with is as follows (depth is in cm):

 Plot()

 {

  P_Sequence("k=0.5",0.5,2,U(-2,2))

  {

   Boundary();

   R_Date("Radio3",924.9619260639960, 22.3517320113017)

   {

    z=115;

   };

   R_Date("Radio2",865.8364138980180, 22.3517320113017)

   {

    z=78;

   };

   R_Date("Radio1",327.7468324823860, 22.3517320113017)

   {

    z=35;

   };

   Date("Top",2015,0)

   {

    z=0;

   };

   Boundary();

  };

 };

As a side-note, I've tried breaking up this profile into smaller sequences/phases, however, the results are not that good because of the sparse dates.

Thanks so much for the help,

Matthew

[Matthew Bolton]

As an update: I did some sensitivity analysis, changing my base k, attempting to get the log(k/k0) mean to zero. Fitting an exponential curve to my points revealed a good estimate of the "realistic" base k. This has been the best (read: only) way I could come up with to determine the optimum k value. I, however, haven't seen anybody else doing this process. I'm probably missing something...

As it turns out, the modelled k for this profile looks something like 0.0685954. Checking the variable k model against a fixed-k model with this value returns comparable (essentially identical) results.

As another side note, extending this event return rate through the profile yields a rather good agreement with the number of defined strata (i.e., large flood events; though the exact "number of events" is dependent on interpretation -- depending on what you define as an "event"). I'm not saying it's the "right" event interval, but this certainly doesn't give me much of a reason to doubt the robustness of OxCal's variable k selection process.

I still look forward to any feedback on this topic -- even if you tell me I'm going about this experiment backwards.

-Matt

[Christopher Ramsey]

Matt

Yes - the value that is modelled in the P_Sequence code is the log10 of k so a -ve value implies a value that is just lower than the one estimated. There is a bottom limit to k because there are events in the sequence which means the number of events per unit m cannot fall below this. If the k value reaches this limit the model is not very different from a simple Sequence. To get the most likely value for k you should be able to take the maximum in the distribution from the P_Sequence and then take this to the power of 10 times the estimated value.

In the default model from the model tool it for cm it sets the k estimate to 1 and there is a U(-2,2) range - so if the peak value in this case was -1 it would imply a k value optimum around 0.1. In many ways it is better to average over results from different appropriate k values so unless there is reason to I would leave the k value variable rather than fixing it. Of course the k value is something you may want to know for other reasons and the model will run faster with fixed k so there may be good reasons to determine it.

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.
> For more options, visit https://groups.google.com/d/optout.

[Matthew Bolton]

Christopher, thanks for your quick reply! It's great to be greeted to the group by the software's originator! 

Anyhow, the value of allowing for variable k is not lost on me. In the meantime, I've been trying to work out exactly what your last message means.

However, I'm still a bit confused with, what I think, should be some simple math. You described "To get the most likely value for k you should be able to take the maximum in the distribution from the P_Sequence and then take this to the power of 10 times the estimated value." For shorthand, I interpret the distribution maximum as the mode of the function, as if it was a histogram. This means (to me) that the modelled k should be as follows: (mode)^10=k. However, this simple conversion doesn't seem to work, particularly since in my case the resultant peak varies wildly with changes in the base k (k0). 

Likewise, for the example where the mode was -1 and the k optimum was 0.01, I still can't recreate these results using any combination of logs and ^10s. Could Prof Ramsey or somebody else perhaps walk me through this simple process in a bit more detail. Attached, are the details for one of my examples where the U(-2,2) range was applied to the base k; in text they are:

k0	Mean log(k) (from OxCal plot)	Modelled Optimum	Sensitivity-predicted k
2.5	-1.56379	?	0.06296

Sorry for the incessant simple questions. Though, it's clear I'm missing a key part of this process. 

Thanks again,

Matt