search

ISAAKiel / quantAAR

R Package - Quantitative Analysis in Archaeology
13 stars 1 forks source link

14C calibration #2

Closed MartinHinz closed 5 years ago

MartinHinz commented 8 years ago

Should 14C calibration be included at all?

If so, should we depend on eg. bchron, or shall we code a calibration by ourselves?

If so, including Bayesian Calibration? Maybe as future milestone?

nevrome commented 8 years ago

Bchron is already on a very good way - nevertheless 14C-calibration and modelling is crucial. Maybe we can set up something in the future.

nevrome commented 8 years ago

I realized that I need a very fast tool to do the calibration for neolithicR. My current solution with Bchron is much to slow.

The last time we talked about this, you described a pretty simple algorithm to perform an equally simple calibration. Do you think we could implement something like this as a Rcpp script quickly?

MartinHinz commented 8 years ago

Sure. I am not sure if Rcpp is really necessary?

A short r-markdown piece below:

Quick Calibration

remark: step one and two are only necessary once, the calibration matrix could be stored

First: Get the cal curve, eg. from the internet (http://www.radiocarbon.org/IntCal13%20files/intcal13.14c)

con <- url("http://www.radiocarbon.org/IntCal13%20files/intcal13.14c")
calcurve<-read.csv(con,skip=11, header=F)
colnames(calcurve)<-c("CAL BP", "14C age","Error","Delta 14C","Sigma")

Second: make a matrix out of it

cal_bp <- calcurve[,1]
cal_matrix <- apply(calcurve,1,function(x){dnorm(cal_bp,mean = x[2], sd = x[3])})

Third: get your 14C date as probability vector

my_date <- c(4000,20)
my_date_prob<-dnorm(cal_bp,mean = my_date[1], sd = my_date[2])

Fourth: Multiply Vector by Matrix

res <- my_date_prob %*% cal_matrix

Fifth: Get the result in more meaningful shape

res.df <- cbind(cal_bp,res[1,])
res.df <- res.df[res.df[,2]>1e-05,]

Sixth: Plot

plot(res.df,type="l")

Seventh: Compare with Bcron:

library(Bchron)

res.bchron<-Bchron::BchronCalibrate(my_date[1],my_date[2],calCurves = "intcal13")
par(mfrow=c(2,1))
plot(res.df,type="l")
plot(res.bchron)
par(mfrow=c(1,1))
MartinHinz commented 8 years ago

quickcal.pdf

nevrome commented 8 years ago

Awesome - thanks! We should really consider to create an own calibration module. The simplicity of this code is impressive and list calibration is probably still the most commonly required task.

MartinHinz commented 5 years ago

In the light of rcarbon and oxcaar, I think there is currently no need for an additional calibration function. So I will close this issue right now.